From 6d8b947d02757278ace320f5d4028b4df28187c5 Mon Sep 17 00:00:00 2001 From: lunamorrow Date: Wed, 27 Mar 2024 17:42:19 +1000 Subject: [PATCH 01/10] changed hbonds --- .../analysis/hydrogenbonds/hbond_analysis.py | 141 ++++++++++-------- 1 file changed, 79 insertions(+), 62 deletions(-) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index e0597e2c3d2..08a4d56e01b 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -406,9 +406,9 @@ def guess_hydrogens(self, Returns ------- - potential_hydrogens: str - String containing the :attr:`resname` and :attr:`name` of all - hydrogen atoms potentially capable of forming hydrogen bonds. + potential_hydrogens: AtomGroup + AtomGroup corresponding to all hydrogen atoms potentially capable + of forming hydrogen bonds. Notes ----- @@ -425,9 +425,9 @@ def guess_hydrogens(self, the selection. Alternatively, this function may be used to quickly generate a - :class:`str` of potential hydrogen atoms involved in hydrogen bonding. - This str may then be modified before being used to set the attribute - :attr:`hydrogens_sel`. + :class:`AtomGroup` of potential hydrogen atoms involved in hydrogen + bonding. This :class:`AtomGroup` may then be modified before being used + to set the attribute :attr:`hydrogens_sel`. .. versionchanged:: 2.4.0 @@ -447,7 +447,7 @@ def guess_hydrogens(self, )) ] - return self._group_categories(hydrogens_ag) + return hydrogens_ag def guess_donors(self, select='all', max_charge=-0.5): """Guesses which atoms could be considered donors in the analysis. Only @@ -465,9 +465,9 @@ def guess_donors(self, select='all', max_charge=-0.5): Returns ------- - potential_donors: str - String containing the :attr:`resname` and :attr:`name` of all atoms - that are potentially capable of forming hydrogen bonds. + potential_donors: AtomGroup + AtomGroup corresponding to all atoms potentially capable of forming + hydrogen bonds. Notes ----- @@ -485,9 +485,9 @@ def guess_donors(self, select='all', max_charge=-0.5): selection. Alternatively, this function may be used to quickly generate a - :class:`str` of potential donor atoms involved in hydrogen bonding. - This :class:`str` may then be modified before being used to set the - attribute :attr:`donors_sel`. + :class:`AtomGroup` of potential donor atoms involved in hydrogen + bonding. This :class:`AtomGroup` may then be modified before being used + to set the attribute :attr:`donors_sel`. .. versionchanged:: 2.4.0 @@ -498,11 +498,11 @@ def guess_donors(self, select='all', max_charge=-0.5): # We need to know `hydrogens_sel` before we can find donors # Use a new variable `hydrogens_sel` so that we do not set # `self.hydrogens_sel` if it is currently `None` + hydrogens_ag = self.guess_hydrogens() if self.hydrogens_sel is None: - hydrogens_sel = self.guess_hydrogens() + hydrogens_sel = self._group_categories(hydrogens_ag) else: hydrogens_sel = self.hydrogens_sel - hydrogens_ag = self.u.select_atoms(hydrogens_sel) # We're using u._topology.bonds rather than u.bonds as it is a million # times faster to access. This is because u.bonds also calculates @@ -521,9 +521,7 @@ def guess_donors(self, select='all', max_charge=-0.5): ) ) - donors_ag = donors_ag[donors_ag.charges < max_charge] - - return self._group_categories(donors_ag) + return donors_ag[donors_ag.charges < max_charge] def guess_acceptors(self, select='all', max_charge=-0.5): """Guesses which atoms could be considered acceptors in the analysis. @@ -542,9 +540,9 @@ def guess_acceptors(self, select='all', max_charge=-0.5): Returns ------- - potential_acceptors: str - String containing the :attr:`resname` and :attr:`name` of all atoms - that potentially capable of forming hydrogen bonds. + potential_acceptors: AtomGroup + AtomGroup corresponding to all atoms potentially capable of forming + hydrogen bonds. Notes ----- @@ -560,9 +558,9 @@ def guess_acceptors(self, select='all', max_charge=-0.5): the selection. Alternatively, this function may be used to quickly generate a - :class:`str` of potential acceptor atoms involved in hydrogen bonding. - This :class:`str` may then be modified before being used to set the - attribute :attr:`acceptors_sel`. + :class:`AtomGroup` of potential acceptor atoms involved in hydrogen + bonding. This :class:`AtomGroup` may then be modified before being used + to set the attribute :attr:`acceptors_sel`. .. versionchanged:: 2.4.0 @@ -570,10 +568,9 @@ def guess_acceptors(self, select='all', max_charge=-0.5): """ - ag = self.u.select_atoms(select) - acceptors_ag = ag[ag.charges < max_charge] + ag = self.u.select_atoms(select, updating=self.update_selections) - return self._group_categories(acceptors_ag) + return ag[ag.charges < max_charge] @staticmethod def _group_categories(group): @@ -621,37 +618,50 @@ def _get_dh_pairs(self): produce a list of donor-hydrogen pairs. """ - # If donors_sel is not provided, use topology to find d-h pairs - if self.donors_sel is None: - - # We're using u._topology.bonds rather than u.bonds as it is a million times faster to access. - # This is because u.bonds also calculates properties of each bond (e.g bond length). - # See https://github.com/MDAnalysis/mdanalysis/issues/2396#issuecomment-596251787 - if not (hasattr(self.u._topology, 'bonds') and len(self.u._topology.bonds.values) != 0): - raise NoDataError('Cannot assign donor-hydrogen pairs via topology as no bond information is present. ' - 'Please either: load a topology file with bond information; use the guess_bonds() ' - 'topology guesser; or set HydrogenBondAnalysis.donors_sel so that a distance cutoff ' - 'can be used.') - - hydrogens = self.u.select_atoms(self.hydrogens_sel) - donors = sum(h.bonded_atoms[0] for h in hydrogens) if hydrogens \ - else AtomGroup([], self.u) - - # Otherwise, use d_h_cutoff as a cutoff distance - else: - - hydrogens = self.u.select_atoms(self.hydrogens_sel) - donors = self.u.select_atoms(self.donors_sel) - donors_indices, hydrogen_indices = capped_distance( - donors.positions, - hydrogens.positions, - max_cutoff=self.d_h_cutoff, - box=self.u.dimensions, - return_distances=False - ).T - - donors = donors[donors_indices] - hydrogens = hydrogens[hydrogen_indices] + # # If donors_sel is not provided, use topology to find d-h pairs + # if self.donors_sel is None: + + # # We're using u._topology.bonds rather than u.bonds as it is a million times faster to access. + # # This is because u.bonds also calculates properties of each bond (e.g bond length). + # # See https://github.com/MDAnalysis/mdanalysis/issues/2396#issuecomment-596251787 + # if not (hasattr(self.u._topology, 'bonds') and len(self.u._topology.bonds.values) != 0): + # raise NoDataError('Cannot assign donor-hydrogen pairs via topology as no bond information is present. ' + # 'Please either: load a topology file with bond information; use the guess_bonds() ' + # 'topology guesser; or set HydrogenBondAnalysis.donors_sel so that a distance cutoff ' + # 'can be used.') + + # hydrogens = self.u.select_atoms(self.hydrogens_sel) + # donors = sum(h.bonded_atoms[0] for h in hydrogens) if hydrogens \ + # else AtomGroup([], self.u) + + # # Otherwise, use d_h_cutoff as a cutoff distance + # else: + + # hydrogens = self.guess_hydrogens() + # donors = self.guess_donors() + # donors_indices, hydrogen_indices = capped_distance( + # donors.positions, + # hydrogens.positions, + # max_cutoff=self.d_h_cutoff, + # box=self.u.dimensions, + # return_distances=False + # ).T + + # donors = donors[donors_indices] + # hydrogens = hydrogens[hydrogen_indices] + + hydrogens = self.hydrogens_ag #NOTE: do I need this? + donors = self.donors_ag #NOTE: do I need this? + donors_indices, hydrogen_indices = capped_distance( + donors.positions, + hydrogens.positions, + max_cutoff=self.d_h_cutoff, + box=self.u.dimensions, + return_distances=False + ).T + + donors = donors[donors_indices] + hydrogens = hydrogens[hydrogen_indices] return donors, hydrogens @@ -699,15 +709,22 @@ def _filter_atoms(self, donors, acceptors): def _prepare(self): self.results.hbonds = [[], [], [], [], [], []] + # Set unpaired acceptor and hydrogen AtomGroups + self.acceptors_ag = self.guess_acceptors() #do i need this? + self.hydrogens_ag = self.guess_hydrogens() + self.donors_ag = self.guess_donors() + # Set atom selections if they have not been provided if self.acceptors_sel is None: - self.acceptors_sel = self.guess_acceptors() + self.acceptors_sel = self._group_categories(self.acceptors_ag) #NOTE: do I need this? if self.hydrogens_sel is None: - self.hydrogens_sel = self.guess_hydrogens() + self.hydrogens_sel = self._group_categories(self.hydrogens_ag) + #NOTE: add in self.donors_sel? Is this not done for a reason? + if self.donors_sel is None: + self.donors_sel = self._group_categories(self.donors_ag) # Select atom groups - self._acceptors = self.u.select_atoms(self.acceptors_sel, - updating=self.update_selections) + self._acceptors = self.guess_acceptors() self._donors, self._hydrogens = self._get_dh_pairs() def _single_frame(self): From 4daea21af12c65e808a60328b0164fb5b3372827 Mon Sep 17 00:00:00 2001 From: lunamorrow Date: Wed, 27 Mar 2024 18:11:00 +1000 Subject: [PATCH 02/10] gsoc-1 --- package/CHANGELOG | 21 +++++++++++++++++++ .../analysis/hydrogenbonds/hbond_analysis.py | 2 +- 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 094ec0091ca..9384caff0bf 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -14,6 +14,27 @@ The rules for this file: ------------------------------------------------------------------------------- +03/27/23 lunamorrow + + * 2.8.1 + +Fixes + +Enhancements + * modernise `HydrogenBondAnalysis` to use `AtomGroups` as objects internally + instead of selection strings + +Changes + * `guess_acceptors`, `guess_hydrogens` and `guess_donators` no longer return + selection strings (return `AtomGroups`instead) + * `_prepare` and `_single_frame` no longer pass selection strings to other + class methods (pass `AtomGroups` instead) + +Deprecations + * Selection string return of `guess_acceptors`, `guess_hydrogens` and + `guess_donators` has been deprecated in favour of return of `AtomGroups` + and will removed in removed in version 3.0.0 (PR #4533) + ??/??/?? IAlibay, HeetVekariya, marinegor, lilyminium, RMeli, ljwoods2, aditya292002, pstaerk, PicoCentauri, BFedder, tyler.je.reddy diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index 08a4d56e01b..c33ae4a7e17 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -710,7 +710,7 @@ def _prepare(self): self.results.hbonds = [[], [], [], [], [], []] # Set unpaired acceptor and hydrogen AtomGroups - self.acceptors_ag = self.guess_acceptors() #do i need this? + self.acceptors_ag = self.guess_acceptors() #NOTE: do I need this? self.hydrogens_ag = self.guess_hydrogens() self.donors_ag = self.guess_donors() From 670781ef85d56e4a6c7f22827b66b02a2c5a3f69 Mon Sep 17 00:00:00 2001 From: lunamorrow Date: Wed, 27 Mar 2024 19:40:03 +1000 Subject: [PATCH 03/10] gsoc-1 --- package/AUTHORS | 1 + .../analysis/hydrogenbonds/hbond_analysis.py | 99 ++++++++++--------- 2 files changed, 53 insertions(+), 47 deletions(-) diff --git a/package/AUTHORS b/package/AUTHORS index 35e4029f776..708226222eb 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -236,6 +236,7 @@ Chronological list of authors 2024 - Aditya Keshari - Philipp Stärk + - Luna Morrow External code ------------- diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index c33ae4a7e17..ca52248e881 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -618,50 +618,52 @@ def _get_dh_pairs(self): produce a list of donor-hydrogen pairs. """ - # # If donors_sel is not provided, use topology to find d-h pairs - # if self.donors_sel is None: - - # # We're using u._topology.bonds rather than u.bonds as it is a million times faster to access. - # # This is because u.bonds also calculates properties of each bond (e.g bond length). - # # See https://github.com/MDAnalysis/mdanalysis/issues/2396#issuecomment-596251787 - # if not (hasattr(self.u._topology, 'bonds') and len(self.u._topology.bonds.values) != 0): - # raise NoDataError('Cannot assign donor-hydrogen pairs via topology as no bond information is present. ' - # 'Please either: load a topology file with bond information; use the guess_bonds() ' - # 'topology guesser; or set HydrogenBondAnalysis.donors_sel so that a distance cutoff ' - # 'can be used.') - - # hydrogens = self.u.select_atoms(self.hydrogens_sel) - # donors = sum(h.bonded_atoms[0] for h in hydrogens) if hydrogens \ - # else AtomGroup([], self.u) - - # # Otherwise, use d_h_cutoff as a cutoff distance - # else: - - # hydrogens = self.guess_hydrogens() - # donors = self.guess_donors() - # donors_indices, hydrogen_indices = capped_distance( - # donors.positions, - # hydrogens.positions, - # max_cutoff=self.d_h_cutoff, - # box=self.u.dimensions, - # return_distances=False - # ).T - - # donors = donors[donors_indices] - # hydrogens = hydrogens[hydrogen_indices] - - hydrogens = self.hydrogens_ag #NOTE: do I need this? - donors = self.donors_ag #NOTE: do I need this? - donors_indices, hydrogen_indices = capped_distance( - donors.positions, - hydrogens.positions, - max_cutoff=self.d_h_cutoff, - box=self.u.dimensions, - return_distances=False - ).T - - donors = donors[donors_indices] - hydrogens = hydrogens[hydrogen_indices] + # If donors_sel is not provided, use topology to find d-h pairs + if self.donors_sel is None: + + # We're using u._topology.bonds rather than u.bonds as it is a million times faster to access. + # This is because u.bonds also calculates properties of each bond (e.g bond length). + # See https://github.com/MDAnalysis/mdanalysis/issues/2396#issuecomment-596251787 + if not (hasattr(self.u._topology, 'bonds') and len(self.u._topology.bonds.values) != 0): + raise NoDataError('Cannot assign donor-hydrogen pairs via topology as no bond information is present. ' + 'Please either: load a topology file with bond information; use the guess_bonds() ' + 'topology guesser; or set HydrogenBondAnalysis.donors_sel so that a distance cutoff ' + 'can be used.') + + hydrogens = self.hydrogens_ag + donors = sum(h.bonded_atoms[0] for h in hydrogens) if hydrogens \ + else AtomGroup([], self.u) + + # Otherwise, use d_h_cutoff as a cutoff distance + else: + + # hydrogens = self.guess_hydrogens() + # donors = self.guess_donors() + hydrogens = self.hydrogens_ag #NOTE: do I need this? + donors = self.donors_ag #NOTE: do I need this? + donors_indices, hydrogen_indices = capped_distance( + donors.positions, + hydrogens.positions, + max_cutoff=self.d_h_cutoff, + box=self.u.dimensions, + return_distances=False + ).T + + donors = donors[donors_indices] + hydrogens = hydrogens[hydrogen_indices] + + # hydrogens = self.hydrogens_ag #NOTE: do I need this? + # donors = self.donors_ag #NOTE: do I need this? + # donors_indices, hydrogen_indices = capped_distance( + # donors.positions, + # hydrogens.positions, + # max_cutoff=self.d_h_cutoff, + # box=self.u.dimensions, + # return_distances=False + # ).T + + # donors = donors[donors_indices] + # hydrogens = hydrogens[hydrogen_indices] return donors, hydrogens @@ -710,17 +712,20 @@ def _prepare(self): self.results.hbonds = [[], [], [], [], [], []] # Set unpaired acceptor and hydrogen AtomGroups - self.acceptors_ag = self.guess_acceptors() #NOTE: do I need this? - self.hydrogens_ag = self.guess_hydrogens() - self.donors_ag = self.guess_donors() + # self.acceptors_ag = self.guess_acceptors() #NOTE: do I need this? + # self.hydrogens_ag = self.guess_hydrogens() + # self.donors_ag = self.guess_donors() # Set atom selections if they have not been provided if self.acceptors_sel is None: + self.acceptors_ag = self.guess_acceptors() self.acceptors_sel = self._group_categories(self.acceptors_ag) #NOTE: do I need this? if self.hydrogens_sel is None: + self.hydrogens_ag = self.guess_hydrogens() self.hydrogens_sel = self._group_categories(self.hydrogens_ag) #NOTE: add in self.donors_sel? Is this not done for a reason? if self.donors_sel is None: + self.donors_ag = self.guess_donors() self.donors_sel = self._group_categories(self.donors_ag) # Select atom groups From 24b5371847c30ec6ad58524fe4f0075baeed81ad Mon Sep 17 00:00:00 2001 From: lunamorrow Date: Tue, 16 Apr 2024 13:10:25 +1000 Subject: [PATCH 04/10] updated-hbond-analysis-with-atomgroups --- .../analysis/hydrogenbonds/hbond_analysis.py | 135 +++++++++--------- .../analysis/test_hydrogenbonds_analysis.py | 10 +- 2 files changed, 77 insertions(+), 68 deletions(-) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index ca52248e881..80d273a61d4 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -406,9 +406,9 @@ def guess_hydrogens(self, Returns ------- - potential_hydrogens: AtomGroup - AtomGroup corresponding to all hydrogen atoms potentially capable - of forming hydrogen bonds. + potential_hydrogens: str + String containing the :attr:`resname` and :attr:`name` of all + hydrogen atoms potentially capable of forming hydrogen bonds. Notes ----- @@ -425,9 +425,9 @@ def guess_hydrogens(self, the selection. Alternatively, this function may be used to quickly generate a - :class:`AtomGroup` of potential hydrogen atoms involved in hydrogen - bonding. This :class:`AtomGroup` may then be modified before being used - to set the attribute :attr:`hydrogens_sel`. + :class:`str` of potential hydrogen atoms involved in hydrogen bonding. + This str may then be modified before being used to set the attribute + :attr:`hydrogens_sel`. .. versionchanged:: 2.4.0 @@ -439,13 +439,17 @@ def guess_hydrogens(self, raise ValueError("min_mass is higher than (or equal to) max_mass") ag = self.u.select_atoms(select) - hydrogens_ag = ag[ - np.logical_and.reduce(( - ag.masses < max_mass, - ag.charges > min_charge, - ag.masses > min_mass, - )) - ] + # hydrogens_ag = ag[ + # np.logical_and.reduce(( + # ag.masses < max_mass, + # ag.charges > min_charge, + # ag.masses > min_mass, + # )) + # ] + hydrogens_ag = ag.select_atoms(f"prop charge > {min_charge} and prop \ + mass > {min_mass} and prop mass < \ + {max_mass}", + updating=self.update_selections) return hydrogens_ag @@ -465,9 +469,9 @@ def guess_donors(self, select='all', max_charge=-0.5): Returns ------- - potential_donors: AtomGroup - AtomGroup corresponding to all atoms potentially capable of forming - hydrogen bonds. + potential_donors: str + String containing the :attr:`resname` and :attr:`name` of all atoms + that are potentially capable of forming hydrogen bonds. Notes ----- @@ -485,9 +489,9 @@ def guess_donors(self, select='all', max_charge=-0.5): selection. Alternatively, this function may be used to quickly generate a - :class:`AtomGroup` of potential donor atoms involved in hydrogen - bonding. This :class:`AtomGroup` may then be modified before being used - to set the attribute :attr:`donors_sel`. + :class:`str` of potential donor atoms involved in hydrogen bonding. + This :class:`str` may then be modified before being used to set the + attribute :attr:`donors_sel`. .. versionchanged:: 2.4.0 @@ -498,11 +502,12 @@ def guess_donors(self, select='all', max_charge=-0.5): # We need to know `hydrogens_sel` before we can find donors # Use a new variable `hydrogens_sel` so that we do not set # `self.hydrogens_sel` if it is currently `None` - hydrogens_ag = self.guess_hydrogens() if self.hydrogens_sel is None: - hydrogens_sel = self._group_categories(hydrogens_ag) + hydrogens_sel = self._group_categories(self.guess_hydrogens()) else: hydrogens_sel = self.hydrogens_sel + hydrogens_ag = self.u.select_atoms(hydrogens_sel, + updating=self.update_selections) # We're using u._topology.bonds rather than u.bonds as it is a million # times faster to access. This is because u.bonds also calculates @@ -521,7 +526,12 @@ def guess_donors(self, select='all', max_charge=-0.5): ) ) - return donors_ag[donors_ag.charges < max_charge] + #donors_ag = donors_ag[donors_ag.charges < max_charge] + + donors_ag = donors_ag.select_atoms(f"prop charge < {max_charge}", + updating=self.update_selections) + + return donors_ag def guess_acceptors(self, select='all', max_charge=-0.5): """Guesses which atoms could be considered acceptors in the analysis. @@ -540,9 +550,9 @@ def guess_acceptors(self, select='all', max_charge=-0.5): Returns ------- - potential_acceptors: AtomGroup - AtomGroup corresponding to all atoms potentially capable of forming - hydrogen bonds. + potential_acceptors: str + String containing the :attr:`resname` and :attr:`name` of all atoms + that potentially capable of forming hydrogen bonds. Notes ----- @@ -558,9 +568,9 @@ def guess_acceptors(self, select='all', max_charge=-0.5): the selection. Alternatively, this function may be used to quickly generate a - :class:`AtomGroup` of potential acceptor atoms involved in hydrogen - bonding. This :class:`AtomGroup` may then be modified before being used - to set the attribute :attr:`acceptors_sel`. + :class:`str` of potential acceptor atoms involved in hydrogen bonding. + This :class:`str` may then be modified before being used to set the + attribute :attr:`acceptors_sel`. .. versionchanged:: 2.4.0 @@ -568,9 +578,12 @@ def guess_acceptors(self, select='all', max_charge=-0.5): """ - ag = self.u.select_atoms(select, updating=self.update_selections) + ag = self.u.select_atoms(select) + #acceptors_ag = ag[ag.charges < max_charge] + acceptors_ag = ag.select_atoms(f"prop charge < {max_charge}", + updating=self.update_selections) - return ag[ag.charges < max_charge] + return acceptors_ag @staticmethod def _group_categories(group): @@ -629,18 +642,24 @@ def _get_dh_pairs(self): 'Please either: load a topology file with bond information; use the guess_bonds() ' 'topology guesser; or set HydrogenBondAnalysis.donors_sel so that a distance cutoff ' 'can be used.') - - hydrogens = self.hydrogens_ag + if self.hydrogens_sel is None: + hydrogens = self.guess_hydrogens() + else: + hydrogens = self.u.select_atoms(self.hydrogens_sel, + updating=self.update_selections) donors = sum(h.bonded_atoms[0] for h in hydrogens) if hydrogens \ else AtomGroup([], self.u) # Otherwise, use d_h_cutoff as a cutoff distance else: - + if self.hydrogens_sel is None: + hydrogens = self.guess_hydrogens() + else: + hydrogens = self.u.select_atoms(self.hydrogens_sel, + updating=self.update_selections) + donors = self.u.select_atoms(self.donors_sel) # hydrogens = self.guess_hydrogens() # donors = self.guess_donors() - hydrogens = self.hydrogens_ag #NOTE: do I need this? - donors = self.donors_ag #NOTE: do I need this? donors_indices, hydrogen_indices = capped_distance( donors.positions, hydrogens.positions, @@ -652,19 +671,6 @@ def _get_dh_pairs(self): donors = donors[donors_indices] hydrogens = hydrogens[hydrogen_indices] - # hydrogens = self.hydrogens_ag #NOTE: do I need this? - # donors = self.donors_ag #NOTE: do I need this? - # donors_indices, hydrogen_indices = capped_distance( - # donors.positions, - # hydrogens.positions, - # max_cutoff=self.d_h_cutoff, - # box=self.u.dimensions, - # return_distances=False - # ).T - - # donors = donors[donors_indices] - # hydrogens = hydrogens[hydrogen_indices] - return donors, hydrogens def _filter_atoms(self, donors, acceptors): @@ -711,25 +717,18 @@ def _filter_atoms(self, donors, acceptors): def _prepare(self): self.results.hbonds = [[], [], [], [], [], []] - # Set unpaired acceptor and hydrogen AtomGroups - # self.acceptors_ag = self.guess_acceptors() #NOTE: do I need this? - # self.hydrogens_ag = self.guess_hydrogens() - # self.donors_ag = self.guess_donors() - # Set atom selections if they have not been provided - if self.acceptors_sel is None: - self.acceptors_ag = self.guess_acceptors() - self.acceptors_sel = self._group_categories(self.acceptors_ag) #NOTE: do I need this? - if self.hydrogens_sel is None: - self.hydrogens_ag = self.guess_hydrogens() - self.hydrogens_sel = self._group_categories(self.hydrogens_ag) - #NOTE: add in self.donors_sel? Is this not done for a reason? - if self.donors_sel is None: - self.donors_ag = self.guess_donors() - self.donors_sel = self._group_categories(self.donors_ag) + # if self.acceptors_sel is None: + # self.acceptors_sel = self.guess_acceptors() + # if self.hydrogens_sel is None: + # self.hydrogens_sel = self.guess_hydrogens() # Select atom groups - self._acceptors = self.guess_acceptors() + if self.acceptors_sel == None: + self._acceptors = self.guess_acceptors() + else: + self._acceptors = self.u.select_atoms(self.acceptors_sel, + updating=self.update_selections) self._donors, self._hydrogens = self._get_dh_pairs() def _single_frame(self): @@ -742,6 +741,8 @@ def _single_frame(self): # find D and A within cutoff distance of one another # min_cutoff = 1.0 as an atom cannot form a hydrogen bond with itself + print(f"D POS: {self._donors.positions}") + print(f"A POS: {self._acceptors.positions}") d_a_indices, d_a_distances = capped_distance( self._donors.positions, self._acceptors.positions, @@ -760,6 +761,9 @@ def _single_frame(self): # Remove D-A pairs more than d_a_cutoff away from one another tmp_donors = self._donors[d_a_indices.T[0]] + print(f"TMP DONORS INITIAL: {tmp_donors}") + print(f"DONORS: {self._donors}") + print(f"DA INDICIES: {d_a_indices.T[0]}") tmp_hydrogens = self._hydrogens[d_a_indices.T[0]] tmp_acceptors = self._acceptors[d_a_indices.T[1]] @@ -792,6 +796,9 @@ def _single_frame(self): # Retrieve atoms, distances and angles of hydrogen bonds hbond_donors = tmp_donors[hbond_indices] + print(f"HBOND DONORS: {hbond_donors}") + print(f"TMP DONORS: {tmp_donors}") + print(f"HBOND INDICIES: {hbond_indices}") hbond_hydrogens = tmp_hydrogens[hbond_indices] hbond_acceptors = tmp_acceptors[hbond_indices] hbond_distances = d_a_distances[hbond_indices] diff --git a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py index b8ea644fc4a..bf9d6d2ef19 100644 --- a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py @@ -227,7 +227,7 @@ def test_no_bond_donor_sel(self, universe): u.add_TopologyAttr('mass', [15.999, 1.008, 1.008] * n_residues) u.add_TopologyAttr('charge', [-1.04, 0.52, 0.52] * n_residues) h = HydrogenBondAnalysis(u, **kwargs) - donors = u.select_atoms(h.guess_donors()) + donors = h.guess_donors() def test_first_hbond(self, hydrogen_bonds): assert len(hydrogen_bonds.results.hbonds) == 2 @@ -554,7 +554,8 @@ def test_guess_donors(self, h): ref_donors = "(resname TIP3 and name OH2)" donors = h.guess_donors(select='all', max_charge=-0.5) - assert donors == ref_donors + donors_sel = h._group_categories(donors) + assert donors_sel == ref_donors class TestHydrogenBondAnalysisTIP3P_GuessHydrogens_NoTopology(object): @@ -586,7 +587,8 @@ def test_guess_hydrogens(self, h): ref_hydrogens = "(resname TIP3 and name H1) or (resname TIP3 and name H2)" hydrogens = h.guess_hydrogens(select='all') - assert hydrogens == ref_hydrogens + hydrogens_sel = h._group_categories(hydrogens) + assert hydrogens_sel == ref_hydrogens pytest.mark.parametrize( "min_mass, max_mass, min_charge", @@ -599,7 +601,7 @@ def test_guess_hydrogens(self, h): def test_guess_hydrogens_empty_selection(self, h): hydrogens = h.guess_hydrogens(select='all', min_charge=1.0) - assert hydrogens == "" + assert len(hydrogens) == 0 def test_guess_hydrogens_min_max_mass(self, h): From 29aa86e6980065417dc1657e6e064743b51fa044 Mon Sep 17 00:00:00 2001 From: lunamorrow Date: Thu, 25 Apr 2024 15:15:51 +1000 Subject: [PATCH 05/10] revised changes, still needs some adjustments and docstring example updates --- package/CHANGELOG | 2 +- .../analysis/hydrogenbonds/hbond_analysis.py | 76 +++++-------------- 2 files changed, 21 insertions(+), 57 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 9384caff0bf..e029083f474 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -14,7 +14,7 @@ The rules for this file: ------------------------------------------------------------------------------- -03/27/23 lunamorrow +04/25/24 lunamorrow * 2.8.1 diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index 80d273a61d4..4118132ad9d 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -108,6 +108,10 @@ Alternatively, :attr:`hydrogens_sel` and :attr:`acceptors_sel` may be generated via the :attr:`guess_hydrogens` and :attr:`guess_acceptors`. This selection strings may then be modified prior to calling :attr:`run`, or a subset of +the universe may be used to guess the atoms. For example, find hydrogens and acceptors belonging to a protein:: + +Alternatively, :attr:`hydrogens_ag` and :attr:`acceptors_ag` may be generated via the :attr:`guess_hydrogens` and +:attr:`guess_acceptors`. This selection strings may then be modified prior to calling :attr:`run`, or a subset of the universe may be used to guess the atoms. For example, find hydrogens and acceptors belonging to a protein:: import MDAnalysis @@ -406,9 +410,9 @@ def guess_hydrogens(self, Returns ------- - potential_hydrogens: str - String containing the :attr:`resname` and :attr:`name` of all - hydrogen atoms potentially capable of forming hydrogen bonds. + potential_hydrogens: AtomGroup + AtomGroup containing all of the hydrogen atoms potentially + capable of forming hydrogen bonds. Notes ----- @@ -424,14 +428,11 @@ def guess_hydrogens(self, If :attr:`hydrogens_sel` is `None`, this function is called to guess the selection. - Alternatively, this function may be used to quickly generate a - :class:`str` of potential hydrogen atoms involved in hydrogen bonding. - This str may then be modified before being used to set the attribute - :attr:`hydrogens_sel`. - .. versionchanged:: 2.4.0 Added ability to use atom types + .. versionchanged:: 3.0.0 + Returns AtomGroup instead of selection string """ @@ -439,13 +440,6 @@ def guess_hydrogens(self, raise ValueError("min_mass is higher than (or equal to) max_mass") ag = self.u.select_atoms(select) - # hydrogens_ag = ag[ - # np.logical_and.reduce(( - # ag.masses < max_mass, - # ag.charges > min_charge, - # ag.masses > min_mass, - # )) - # ] hydrogens_ag = ag.select_atoms(f"prop charge > {min_charge} and prop \ mass > {min_mass} and prop mass < \ {max_mass}", @@ -469,9 +463,9 @@ def guess_donors(self, select='all', max_charge=-0.5): Returns ------- - potential_donors: str - String containing the :attr:`resname` and :attr:`name` of all atoms - that are potentially capable of forming hydrogen bonds. + potential_donors: AtomGroup + AtomGroup containing all of the atoms potentially capable of + forming hydrogen bonds. Notes ----- @@ -488,14 +482,11 @@ def guess_donors(self, select='all', max_charge=-0.5): have bonding information, this function is called to guess the selection. - Alternatively, this function may be used to quickly generate a - :class:`str` of potential donor atoms involved in hydrogen bonding. - This :class:`str` may then be modified before being used to set the - attribute :attr:`donors_sel`. - .. versionchanged:: 2.4.0 Added ability to use atom types + .. versionchanged:: 3.0.0 + Returns AtomGroup instead of selection string """ @@ -526,8 +517,6 @@ def guess_donors(self, select='all', max_charge=-0.5): ) ) - #donors_ag = donors_ag[donors_ag.charges < max_charge] - donors_ag = donors_ag.select_atoms(f"prop charge < {max_charge}", updating=self.update_selections) @@ -550,9 +539,9 @@ def guess_acceptors(self, select='all', max_charge=-0.5): Returns ------- - potential_acceptors: str - String containing the :attr:`resname` and :attr:`name` of all atoms - that potentially capable of forming hydrogen bonds. + potential_acceptors: AtomGroup + AtomGroup containing all of the atoms potentially capable of + forming hydrogen bonds. Notes ----- @@ -567,14 +556,11 @@ def guess_acceptors(self, select='all', max_charge=-0.5): If :attr:`acceptors_sel` is `None`, this function is called to guess the selection. - Alternatively, this function may be used to quickly generate a - :class:`str` of potential acceptor atoms involved in hydrogen bonding. - This :class:`str` may then be modified before being used to set the - attribute :attr:`acceptors_sel`. - .. versionchanged:: 2.4.0 Added ability to use atom types + .. versionchanged:: 3.0.0 + Returns AtomGroup instead of selection string """ @@ -658,8 +644,6 @@ def _get_dh_pairs(self): hydrogens = self.u.select_atoms(self.hydrogens_sel, updating=self.update_selections) donors = self.u.select_atoms(self.donors_sel) - # hydrogens = self.guess_hydrogens() - # donors = self.guess_donors() donors_indices, hydrogen_indices = capped_distance( donors.positions, hydrogens.positions, @@ -717,14 +701,8 @@ def _filter_atoms(self, donors, acceptors): def _prepare(self): self.results.hbonds = [[], [], [], [], [], []] - # Set atom selections if they have not been provided - # if self.acceptors_sel is None: - # self.acceptors_sel = self.guess_acceptors() - # if self.hydrogens_sel is None: - # self.hydrogens_sel = self.guess_hydrogens() - # Select atom groups - if self.acceptors_sel == None: + if self.acceptors_sel is None: self._acceptors = self.guess_acceptors() else: self._acceptors = self.u.select_atoms(self.acceptors_sel, @@ -735,14 +713,6 @@ def _single_frame(self): box = self._ts.dimensions - # Update donor-hydrogen pairs if necessary - if self.update_selections: - self._donors, self._hydrogens = self._get_dh_pairs() - - # find D and A within cutoff distance of one another - # min_cutoff = 1.0 as an atom cannot form a hydrogen bond with itself - print(f"D POS: {self._donors.positions}") - print(f"A POS: {self._acceptors.positions}") d_a_indices, d_a_distances = capped_distance( self._donors.positions, self._acceptors.positions, @@ -761,9 +731,6 @@ def _single_frame(self): # Remove D-A pairs more than d_a_cutoff away from one another tmp_donors = self._donors[d_a_indices.T[0]] - print(f"TMP DONORS INITIAL: {tmp_donors}") - print(f"DONORS: {self._donors}") - print(f"DA INDICIES: {d_a_indices.T[0]}") tmp_hydrogens = self._hydrogens[d_a_indices.T[0]] tmp_acceptors = self._acceptors[d_a_indices.T[1]] @@ -796,9 +763,6 @@ def _single_frame(self): # Retrieve atoms, distances and angles of hydrogen bonds hbond_donors = tmp_donors[hbond_indices] - print(f"HBOND DONORS: {hbond_donors}") - print(f"TMP DONORS: {tmp_donors}") - print(f"HBOND INDICIES: {hbond_indices}") hbond_hydrogens = tmp_hydrogens[hbond_indices] hbond_acceptors = tmp_acceptors[hbond_indices] hbond_distances = d_a_distances[hbond_indices] From 466d6c008a355ca0c7a475e566d1a236ccade142 Mon Sep 17 00:00:00 2001 From: lunamorrow Date: Wed, 27 Mar 2024 17:42:19 +1000 Subject: [PATCH 06/10] changed hbonds --- .../analysis/hydrogenbonds/hbond_analysis.py | 141 ++++++++++-------- 1 file changed, 79 insertions(+), 62 deletions(-) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index e0597e2c3d2..08a4d56e01b 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -406,9 +406,9 @@ def guess_hydrogens(self, Returns ------- - potential_hydrogens: str - String containing the :attr:`resname` and :attr:`name` of all - hydrogen atoms potentially capable of forming hydrogen bonds. + potential_hydrogens: AtomGroup + AtomGroup corresponding to all hydrogen atoms potentially capable + of forming hydrogen bonds. Notes ----- @@ -425,9 +425,9 @@ def guess_hydrogens(self, the selection. Alternatively, this function may be used to quickly generate a - :class:`str` of potential hydrogen atoms involved in hydrogen bonding. - This str may then be modified before being used to set the attribute - :attr:`hydrogens_sel`. + :class:`AtomGroup` of potential hydrogen atoms involved in hydrogen + bonding. This :class:`AtomGroup` may then be modified before being used + to set the attribute :attr:`hydrogens_sel`. .. versionchanged:: 2.4.0 @@ -447,7 +447,7 @@ def guess_hydrogens(self, )) ] - return self._group_categories(hydrogens_ag) + return hydrogens_ag def guess_donors(self, select='all', max_charge=-0.5): """Guesses which atoms could be considered donors in the analysis. Only @@ -465,9 +465,9 @@ def guess_donors(self, select='all', max_charge=-0.5): Returns ------- - potential_donors: str - String containing the :attr:`resname` and :attr:`name` of all atoms - that are potentially capable of forming hydrogen bonds. + potential_donors: AtomGroup + AtomGroup corresponding to all atoms potentially capable of forming + hydrogen bonds. Notes ----- @@ -485,9 +485,9 @@ def guess_donors(self, select='all', max_charge=-0.5): selection. Alternatively, this function may be used to quickly generate a - :class:`str` of potential donor atoms involved in hydrogen bonding. - This :class:`str` may then be modified before being used to set the - attribute :attr:`donors_sel`. + :class:`AtomGroup` of potential donor atoms involved in hydrogen + bonding. This :class:`AtomGroup` may then be modified before being used + to set the attribute :attr:`donors_sel`. .. versionchanged:: 2.4.0 @@ -498,11 +498,11 @@ def guess_donors(self, select='all', max_charge=-0.5): # We need to know `hydrogens_sel` before we can find donors # Use a new variable `hydrogens_sel` so that we do not set # `self.hydrogens_sel` if it is currently `None` + hydrogens_ag = self.guess_hydrogens() if self.hydrogens_sel is None: - hydrogens_sel = self.guess_hydrogens() + hydrogens_sel = self._group_categories(hydrogens_ag) else: hydrogens_sel = self.hydrogens_sel - hydrogens_ag = self.u.select_atoms(hydrogens_sel) # We're using u._topology.bonds rather than u.bonds as it is a million # times faster to access. This is because u.bonds also calculates @@ -521,9 +521,7 @@ def guess_donors(self, select='all', max_charge=-0.5): ) ) - donors_ag = donors_ag[donors_ag.charges < max_charge] - - return self._group_categories(donors_ag) + return donors_ag[donors_ag.charges < max_charge] def guess_acceptors(self, select='all', max_charge=-0.5): """Guesses which atoms could be considered acceptors in the analysis. @@ -542,9 +540,9 @@ def guess_acceptors(self, select='all', max_charge=-0.5): Returns ------- - potential_acceptors: str - String containing the :attr:`resname` and :attr:`name` of all atoms - that potentially capable of forming hydrogen bonds. + potential_acceptors: AtomGroup + AtomGroup corresponding to all atoms potentially capable of forming + hydrogen bonds. Notes ----- @@ -560,9 +558,9 @@ def guess_acceptors(self, select='all', max_charge=-0.5): the selection. Alternatively, this function may be used to quickly generate a - :class:`str` of potential acceptor atoms involved in hydrogen bonding. - This :class:`str` may then be modified before being used to set the - attribute :attr:`acceptors_sel`. + :class:`AtomGroup` of potential acceptor atoms involved in hydrogen + bonding. This :class:`AtomGroup` may then be modified before being used + to set the attribute :attr:`acceptors_sel`. .. versionchanged:: 2.4.0 @@ -570,10 +568,9 @@ def guess_acceptors(self, select='all', max_charge=-0.5): """ - ag = self.u.select_atoms(select) - acceptors_ag = ag[ag.charges < max_charge] + ag = self.u.select_atoms(select, updating=self.update_selections) - return self._group_categories(acceptors_ag) + return ag[ag.charges < max_charge] @staticmethod def _group_categories(group): @@ -621,37 +618,50 @@ def _get_dh_pairs(self): produce a list of donor-hydrogen pairs. """ - # If donors_sel is not provided, use topology to find d-h pairs - if self.donors_sel is None: - - # We're using u._topology.bonds rather than u.bonds as it is a million times faster to access. - # This is because u.bonds also calculates properties of each bond (e.g bond length). - # See https://github.com/MDAnalysis/mdanalysis/issues/2396#issuecomment-596251787 - if not (hasattr(self.u._topology, 'bonds') and len(self.u._topology.bonds.values) != 0): - raise NoDataError('Cannot assign donor-hydrogen pairs via topology as no bond information is present. ' - 'Please either: load a topology file with bond information; use the guess_bonds() ' - 'topology guesser; or set HydrogenBondAnalysis.donors_sel so that a distance cutoff ' - 'can be used.') - - hydrogens = self.u.select_atoms(self.hydrogens_sel) - donors = sum(h.bonded_atoms[0] for h in hydrogens) if hydrogens \ - else AtomGroup([], self.u) - - # Otherwise, use d_h_cutoff as a cutoff distance - else: - - hydrogens = self.u.select_atoms(self.hydrogens_sel) - donors = self.u.select_atoms(self.donors_sel) - donors_indices, hydrogen_indices = capped_distance( - donors.positions, - hydrogens.positions, - max_cutoff=self.d_h_cutoff, - box=self.u.dimensions, - return_distances=False - ).T - - donors = donors[donors_indices] - hydrogens = hydrogens[hydrogen_indices] + # # If donors_sel is not provided, use topology to find d-h pairs + # if self.donors_sel is None: + + # # We're using u._topology.bonds rather than u.bonds as it is a million times faster to access. + # # This is because u.bonds also calculates properties of each bond (e.g bond length). + # # See https://github.com/MDAnalysis/mdanalysis/issues/2396#issuecomment-596251787 + # if not (hasattr(self.u._topology, 'bonds') and len(self.u._topology.bonds.values) != 0): + # raise NoDataError('Cannot assign donor-hydrogen pairs via topology as no bond information is present. ' + # 'Please either: load a topology file with bond information; use the guess_bonds() ' + # 'topology guesser; or set HydrogenBondAnalysis.donors_sel so that a distance cutoff ' + # 'can be used.') + + # hydrogens = self.u.select_atoms(self.hydrogens_sel) + # donors = sum(h.bonded_atoms[0] for h in hydrogens) if hydrogens \ + # else AtomGroup([], self.u) + + # # Otherwise, use d_h_cutoff as a cutoff distance + # else: + + # hydrogens = self.guess_hydrogens() + # donors = self.guess_donors() + # donors_indices, hydrogen_indices = capped_distance( + # donors.positions, + # hydrogens.positions, + # max_cutoff=self.d_h_cutoff, + # box=self.u.dimensions, + # return_distances=False + # ).T + + # donors = donors[donors_indices] + # hydrogens = hydrogens[hydrogen_indices] + + hydrogens = self.hydrogens_ag #NOTE: do I need this? + donors = self.donors_ag #NOTE: do I need this? + donors_indices, hydrogen_indices = capped_distance( + donors.positions, + hydrogens.positions, + max_cutoff=self.d_h_cutoff, + box=self.u.dimensions, + return_distances=False + ).T + + donors = donors[donors_indices] + hydrogens = hydrogens[hydrogen_indices] return donors, hydrogens @@ -699,15 +709,22 @@ def _filter_atoms(self, donors, acceptors): def _prepare(self): self.results.hbonds = [[], [], [], [], [], []] + # Set unpaired acceptor and hydrogen AtomGroups + self.acceptors_ag = self.guess_acceptors() #do i need this? + self.hydrogens_ag = self.guess_hydrogens() + self.donors_ag = self.guess_donors() + # Set atom selections if they have not been provided if self.acceptors_sel is None: - self.acceptors_sel = self.guess_acceptors() + self.acceptors_sel = self._group_categories(self.acceptors_ag) #NOTE: do I need this? if self.hydrogens_sel is None: - self.hydrogens_sel = self.guess_hydrogens() + self.hydrogens_sel = self._group_categories(self.hydrogens_ag) + #NOTE: add in self.donors_sel? Is this not done for a reason? + if self.donors_sel is None: + self.donors_sel = self._group_categories(self.donors_ag) # Select atom groups - self._acceptors = self.u.select_atoms(self.acceptors_sel, - updating=self.update_selections) + self._acceptors = self.guess_acceptors() self._donors, self._hydrogens = self._get_dh_pairs() def _single_frame(self): From 7a514dc1200a2c30d8ccfa3aac6e6a624a1f0671 Mon Sep 17 00:00:00 2001 From: lunamorrow Date: Wed, 27 Mar 2024 18:11:00 +1000 Subject: [PATCH 07/10] gsoc-1 --- package/CHANGELOG | 21 +++++++++++++++++++ .../analysis/hydrogenbonds/hbond_analysis.py | 2 +- 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index e59f1b068d0..44a246745d5 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -14,6 +14,27 @@ The rules for this file: ------------------------------------------------------------------------------- +03/27/23 lunamorrow + + * 2.8.1 + +Fixes + +Enhancements + * modernise `HydrogenBondAnalysis` to use `AtomGroups` as objects internally + instead of selection strings + +Changes + * `guess_acceptors`, `guess_hydrogens` and `guess_donators` no longer return + selection strings (return `AtomGroups`instead) + * `_prepare` and `_single_frame` no longer pass selection strings to other + class methods (pass `AtomGroups` instead) + +Deprecations + * Selection string return of `guess_acceptors`, `guess_hydrogens` and + `guess_donators` has been deprecated in favour of return of `AtomGroups` + and will removed in removed in version 3.0.0 (PR #4533) + ??/??/?? IAlibay, HeetVekariya, marinegor, lilyminium, RMeli, ljwoods2, aditya292002, pstaerk, PicoCentauri, BFedder, tyler.je.reddy, SampurnaM, leonwehrhan, kainszs, orionarcher diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index 08a4d56e01b..c33ae4a7e17 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -710,7 +710,7 @@ def _prepare(self): self.results.hbonds = [[], [], [], [], [], []] # Set unpaired acceptor and hydrogen AtomGroups - self.acceptors_ag = self.guess_acceptors() #do i need this? + self.acceptors_ag = self.guess_acceptors() #NOTE: do I need this? self.hydrogens_ag = self.guess_hydrogens() self.donors_ag = self.guess_donors() From 4f8d1c1a252fac2b6acbd35f9806ccf5d7d9c693 Mon Sep 17 00:00:00 2001 From: lunamorrow Date: Wed, 27 Mar 2024 19:40:03 +1000 Subject: [PATCH 08/10] gsoc-1 --- package/AUTHORS | 3 +- .../analysis/hydrogenbonds/hbond_analysis.py | 99 ++++++++++--------- 2 files changed, 53 insertions(+), 49 deletions(-) diff --git a/package/AUTHORS b/package/AUTHORS index fb5d5689468..b43de645e47 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -240,8 +240,7 @@ Chronological list of authors - Sampurna Mukherjee - Leon Wehrhan - Valerij Talagayev - - + - Luna Morrow External code ------------- diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index c33ae4a7e17..ca52248e881 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -618,50 +618,52 @@ def _get_dh_pairs(self): produce a list of donor-hydrogen pairs. """ - # # If donors_sel is not provided, use topology to find d-h pairs - # if self.donors_sel is None: - - # # We're using u._topology.bonds rather than u.bonds as it is a million times faster to access. - # # This is because u.bonds also calculates properties of each bond (e.g bond length). - # # See https://github.com/MDAnalysis/mdanalysis/issues/2396#issuecomment-596251787 - # if not (hasattr(self.u._topology, 'bonds') and len(self.u._topology.bonds.values) != 0): - # raise NoDataError('Cannot assign donor-hydrogen pairs via topology as no bond information is present. ' - # 'Please either: load a topology file with bond information; use the guess_bonds() ' - # 'topology guesser; or set HydrogenBondAnalysis.donors_sel so that a distance cutoff ' - # 'can be used.') - - # hydrogens = self.u.select_atoms(self.hydrogens_sel) - # donors = sum(h.bonded_atoms[0] for h in hydrogens) if hydrogens \ - # else AtomGroup([], self.u) - - # # Otherwise, use d_h_cutoff as a cutoff distance - # else: - - # hydrogens = self.guess_hydrogens() - # donors = self.guess_donors() - # donors_indices, hydrogen_indices = capped_distance( - # donors.positions, - # hydrogens.positions, - # max_cutoff=self.d_h_cutoff, - # box=self.u.dimensions, - # return_distances=False - # ).T - - # donors = donors[donors_indices] - # hydrogens = hydrogens[hydrogen_indices] - - hydrogens = self.hydrogens_ag #NOTE: do I need this? - donors = self.donors_ag #NOTE: do I need this? - donors_indices, hydrogen_indices = capped_distance( - donors.positions, - hydrogens.positions, - max_cutoff=self.d_h_cutoff, - box=self.u.dimensions, - return_distances=False - ).T - - donors = donors[donors_indices] - hydrogens = hydrogens[hydrogen_indices] + # If donors_sel is not provided, use topology to find d-h pairs + if self.donors_sel is None: + + # We're using u._topology.bonds rather than u.bonds as it is a million times faster to access. + # This is because u.bonds also calculates properties of each bond (e.g bond length). + # See https://github.com/MDAnalysis/mdanalysis/issues/2396#issuecomment-596251787 + if not (hasattr(self.u._topology, 'bonds') and len(self.u._topology.bonds.values) != 0): + raise NoDataError('Cannot assign donor-hydrogen pairs via topology as no bond information is present. ' + 'Please either: load a topology file with bond information; use the guess_bonds() ' + 'topology guesser; or set HydrogenBondAnalysis.donors_sel so that a distance cutoff ' + 'can be used.') + + hydrogens = self.hydrogens_ag + donors = sum(h.bonded_atoms[0] for h in hydrogens) if hydrogens \ + else AtomGroup([], self.u) + + # Otherwise, use d_h_cutoff as a cutoff distance + else: + + # hydrogens = self.guess_hydrogens() + # donors = self.guess_donors() + hydrogens = self.hydrogens_ag #NOTE: do I need this? + donors = self.donors_ag #NOTE: do I need this? + donors_indices, hydrogen_indices = capped_distance( + donors.positions, + hydrogens.positions, + max_cutoff=self.d_h_cutoff, + box=self.u.dimensions, + return_distances=False + ).T + + donors = donors[donors_indices] + hydrogens = hydrogens[hydrogen_indices] + + # hydrogens = self.hydrogens_ag #NOTE: do I need this? + # donors = self.donors_ag #NOTE: do I need this? + # donors_indices, hydrogen_indices = capped_distance( + # donors.positions, + # hydrogens.positions, + # max_cutoff=self.d_h_cutoff, + # box=self.u.dimensions, + # return_distances=False + # ).T + + # donors = donors[donors_indices] + # hydrogens = hydrogens[hydrogen_indices] return donors, hydrogens @@ -710,17 +712,20 @@ def _prepare(self): self.results.hbonds = [[], [], [], [], [], []] # Set unpaired acceptor and hydrogen AtomGroups - self.acceptors_ag = self.guess_acceptors() #NOTE: do I need this? - self.hydrogens_ag = self.guess_hydrogens() - self.donors_ag = self.guess_donors() + # self.acceptors_ag = self.guess_acceptors() #NOTE: do I need this? + # self.hydrogens_ag = self.guess_hydrogens() + # self.donors_ag = self.guess_donors() # Set atom selections if they have not been provided if self.acceptors_sel is None: + self.acceptors_ag = self.guess_acceptors() self.acceptors_sel = self._group_categories(self.acceptors_ag) #NOTE: do I need this? if self.hydrogens_sel is None: + self.hydrogens_ag = self.guess_hydrogens() self.hydrogens_sel = self._group_categories(self.hydrogens_ag) #NOTE: add in self.donors_sel? Is this not done for a reason? if self.donors_sel is None: + self.donors_ag = self.guess_donors() self.donors_sel = self._group_categories(self.donors_ag) # Select atom groups From f0bd521236d25b04e6e9b66ba387b77e5d388ed8 Mon Sep 17 00:00:00 2001 From: lunamorrow Date: Tue, 16 Apr 2024 13:10:25 +1000 Subject: [PATCH 09/10] updated-hbond-analysis-with-atomgroups --- .../analysis/hydrogenbonds/hbond_analysis.py | 135 +++++++++--------- .../analysis/test_hydrogenbonds_analysis.py | 10 +- 2 files changed, 77 insertions(+), 68 deletions(-) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index ca52248e881..80d273a61d4 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -406,9 +406,9 @@ def guess_hydrogens(self, Returns ------- - potential_hydrogens: AtomGroup - AtomGroup corresponding to all hydrogen atoms potentially capable - of forming hydrogen bonds. + potential_hydrogens: str + String containing the :attr:`resname` and :attr:`name` of all + hydrogen atoms potentially capable of forming hydrogen bonds. Notes ----- @@ -425,9 +425,9 @@ def guess_hydrogens(self, the selection. Alternatively, this function may be used to quickly generate a - :class:`AtomGroup` of potential hydrogen atoms involved in hydrogen - bonding. This :class:`AtomGroup` may then be modified before being used - to set the attribute :attr:`hydrogens_sel`. + :class:`str` of potential hydrogen atoms involved in hydrogen bonding. + This str may then be modified before being used to set the attribute + :attr:`hydrogens_sel`. .. versionchanged:: 2.4.0 @@ -439,13 +439,17 @@ def guess_hydrogens(self, raise ValueError("min_mass is higher than (or equal to) max_mass") ag = self.u.select_atoms(select) - hydrogens_ag = ag[ - np.logical_and.reduce(( - ag.masses < max_mass, - ag.charges > min_charge, - ag.masses > min_mass, - )) - ] + # hydrogens_ag = ag[ + # np.logical_and.reduce(( + # ag.masses < max_mass, + # ag.charges > min_charge, + # ag.masses > min_mass, + # )) + # ] + hydrogens_ag = ag.select_atoms(f"prop charge > {min_charge} and prop \ + mass > {min_mass} and prop mass < \ + {max_mass}", + updating=self.update_selections) return hydrogens_ag @@ -465,9 +469,9 @@ def guess_donors(self, select='all', max_charge=-0.5): Returns ------- - potential_donors: AtomGroup - AtomGroup corresponding to all atoms potentially capable of forming - hydrogen bonds. + potential_donors: str + String containing the :attr:`resname` and :attr:`name` of all atoms + that are potentially capable of forming hydrogen bonds. Notes ----- @@ -485,9 +489,9 @@ def guess_donors(self, select='all', max_charge=-0.5): selection. Alternatively, this function may be used to quickly generate a - :class:`AtomGroup` of potential donor atoms involved in hydrogen - bonding. This :class:`AtomGroup` may then be modified before being used - to set the attribute :attr:`donors_sel`. + :class:`str` of potential donor atoms involved in hydrogen bonding. + This :class:`str` may then be modified before being used to set the + attribute :attr:`donors_sel`. .. versionchanged:: 2.4.0 @@ -498,11 +502,12 @@ def guess_donors(self, select='all', max_charge=-0.5): # We need to know `hydrogens_sel` before we can find donors # Use a new variable `hydrogens_sel` so that we do not set # `self.hydrogens_sel` if it is currently `None` - hydrogens_ag = self.guess_hydrogens() if self.hydrogens_sel is None: - hydrogens_sel = self._group_categories(hydrogens_ag) + hydrogens_sel = self._group_categories(self.guess_hydrogens()) else: hydrogens_sel = self.hydrogens_sel + hydrogens_ag = self.u.select_atoms(hydrogens_sel, + updating=self.update_selections) # We're using u._topology.bonds rather than u.bonds as it is a million # times faster to access. This is because u.bonds also calculates @@ -521,7 +526,12 @@ def guess_donors(self, select='all', max_charge=-0.5): ) ) - return donors_ag[donors_ag.charges < max_charge] + #donors_ag = donors_ag[donors_ag.charges < max_charge] + + donors_ag = donors_ag.select_atoms(f"prop charge < {max_charge}", + updating=self.update_selections) + + return donors_ag def guess_acceptors(self, select='all', max_charge=-0.5): """Guesses which atoms could be considered acceptors in the analysis. @@ -540,9 +550,9 @@ def guess_acceptors(self, select='all', max_charge=-0.5): Returns ------- - potential_acceptors: AtomGroup - AtomGroup corresponding to all atoms potentially capable of forming - hydrogen bonds. + potential_acceptors: str + String containing the :attr:`resname` and :attr:`name` of all atoms + that potentially capable of forming hydrogen bonds. Notes ----- @@ -558,9 +568,9 @@ def guess_acceptors(self, select='all', max_charge=-0.5): the selection. Alternatively, this function may be used to quickly generate a - :class:`AtomGroup` of potential acceptor atoms involved in hydrogen - bonding. This :class:`AtomGroup` may then be modified before being used - to set the attribute :attr:`acceptors_sel`. + :class:`str` of potential acceptor atoms involved in hydrogen bonding. + This :class:`str` may then be modified before being used to set the + attribute :attr:`acceptors_sel`. .. versionchanged:: 2.4.0 @@ -568,9 +578,12 @@ def guess_acceptors(self, select='all', max_charge=-0.5): """ - ag = self.u.select_atoms(select, updating=self.update_selections) + ag = self.u.select_atoms(select) + #acceptors_ag = ag[ag.charges < max_charge] + acceptors_ag = ag.select_atoms(f"prop charge < {max_charge}", + updating=self.update_selections) - return ag[ag.charges < max_charge] + return acceptors_ag @staticmethod def _group_categories(group): @@ -629,18 +642,24 @@ def _get_dh_pairs(self): 'Please either: load a topology file with bond information; use the guess_bonds() ' 'topology guesser; or set HydrogenBondAnalysis.donors_sel so that a distance cutoff ' 'can be used.') - - hydrogens = self.hydrogens_ag + if self.hydrogens_sel is None: + hydrogens = self.guess_hydrogens() + else: + hydrogens = self.u.select_atoms(self.hydrogens_sel, + updating=self.update_selections) donors = sum(h.bonded_atoms[0] for h in hydrogens) if hydrogens \ else AtomGroup([], self.u) # Otherwise, use d_h_cutoff as a cutoff distance else: - + if self.hydrogens_sel is None: + hydrogens = self.guess_hydrogens() + else: + hydrogens = self.u.select_atoms(self.hydrogens_sel, + updating=self.update_selections) + donors = self.u.select_atoms(self.donors_sel) # hydrogens = self.guess_hydrogens() # donors = self.guess_donors() - hydrogens = self.hydrogens_ag #NOTE: do I need this? - donors = self.donors_ag #NOTE: do I need this? donors_indices, hydrogen_indices = capped_distance( donors.positions, hydrogens.positions, @@ -652,19 +671,6 @@ def _get_dh_pairs(self): donors = donors[donors_indices] hydrogens = hydrogens[hydrogen_indices] - # hydrogens = self.hydrogens_ag #NOTE: do I need this? - # donors = self.donors_ag #NOTE: do I need this? - # donors_indices, hydrogen_indices = capped_distance( - # donors.positions, - # hydrogens.positions, - # max_cutoff=self.d_h_cutoff, - # box=self.u.dimensions, - # return_distances=False - # ).T - - # donors = donors[donors_indices] - # hydrogens = hydrogens[hydrogen_indices] - return donors, hydrogens def _filter_atoms(self, donors, acceptors): @@ -711,25 +717,18 @@ def _filter_atoms(self, donors, acceptors): def _prepare(self): self.results.hbonds = [[], [], [], [], [], []] - # Set unpaired acceptor and hydrogen AtomGroups - # self.acceptors_ag = self.guess_acceptors() #NOTE: do I need this? - # self.hydrogens_ag = self.guess_hydrogens() - # self.donors_ag = self.guess_donors() - # Set atom selections if they have not been provided - if self.acceptors_sel is None: - self.acceptors_ag = self.guess_acceptors() - self.acceptors_sel = self._group_categories(self.acceptors_ag) #NOTE: do I need this? - if self.hydrogens_sel is None: - self.hydrogens_ag = self.guess_hydrogens() - self.hydrogens_sel = self._group_categories(self.hydrogens_ag) - #NOTE: add in self.donors_sel? Is this not done for a reason? - if self.donors_sel is None: - self.donors_ag = self.guess_donors() - self.donors_sel = self._group_categories(self.donors_ag) + # if self.acceptors_sel is None: + # self.acceptors_sel = self.guess_acceptors() + # if self.hydrogens_sel is None: + # self.hydrogens_sel = self.guess_hydrogens() # Select atom groups - self._acceptors = self.guess_acceptors() + if self.acceptors_sel == None: + self._acceptors = self.guess_acceptors() + else: + self._acceptors = self.u.select_atoms(self.acceptors_sel, + updating=self.update_selections) self._donors, self._hydrogens = self._get_dh_pairs() def _single_frame(self): @@ -742,6 +741,8 @@ def _single_frame(self): # find D and A within cutoff distance of one another # min_cutoff = 1.0 as an atom cannot form a hydrogen bond with itself + print(f"D POS: {self._donors.positions}") + print(f"A POS: {self._acceptors.positions}") d_a_indices, d_a_distances = capped_distance( self._donors.positions, self._acceptors.positions, @@ -760,6 +761,9 @@ def _single_frame(self): # Remove D-A pairs more than d_a_cutoff away from one another tmp_donors = self._donors[d_a_indices.T[0]] + print(f"TMP DONORS INITIAL: {tmp_donors}") + print(f"DONORS: {self._donors}") + print(f"DA INDICIES: {d_a_indices.T[0]}") tmp_hydrogens = self._hydrogens[d_a_indices.T[0]] tmp_acceptors = self._acceptors[d_a_indices.T[1]] @@ -792,6 +796,9 @@ def _single_frame(self): # Retrieve atoms, distances and angles of hydrogen bonds hbond_donors = tmp_donors[hbond_indices] + print(f"HBOND DONORS: {hbond_donors}") + print(f"TMP DONORS: {tmp_donors}") + print(f"HBOND INDICIES: {hbond_indices}") hbond_hydrogens = tmp_hydrogens[hbond_indices] hbond_acceptors = tmp_acceptors[hbond_indices] hbond_distances = d_a_distances[hbond_indices] diff --git a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py index b8ea644fc4a..bf9d6d2ef19 100644 --- a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py @@ -227,7 +227,7 @@ def test_no_bond_donor_sel(self, universe): u.add_TopologyAttr('mass', [15.999, 1.008, 1.008] * n_residues) u.add_TopologyAttr('charge', [-1.04, 0.52, 0.52] * n_residues) h = HydrogenBondAnalysis(u, **kwargs) - donors = u.select_atoms(h.guess_donors()) + donors = h.guess_donors() def test_first_hbond(self, hydrogen_bonds): assert len(hydrogen_bonds.results.hbonds) == 2 @@ -554,7 +554,8 @@ def test_guess_donors(self, h): ref_donors = "(resname TIP3 and name OH2)" donors = h.guess_donors(select='all', max_charge=-0.5) - assert donors == ref_donors + donors_sel = h._group_categories(donors) + assert donors_sel == ref_donors class TestHydrogenBondAnalysisTIP3P_GuessHydrogens_NoTopology(object): @@ -586,7 +587,8 @@ def test_guess_hydrogens(self, h): ref_hydrogens = "(resname TIP3 and name H1) or (resname TIP3 and name H2)" hydrogens = h.guess_hydrogens(select='all') - assert hydrogens == ref_hydrogens + hydrogens_sel = h._group_categories(hydrogens) + assert hydrogens_sel == ref_hydrogens pytest.mark.parametrize( "min_mass, max_mass, min_charge", @@ -599,7 +601,7 @@ def test_guess_hydrogens(self, h): def test_guess_hydrogens_empty_selection(self, h): hydrogens = h.guess_hydrogens(select='all', min_charge=1.0) - assert hydrogens == "" + assert len(hydrogens) == 0 def test_guess_hydrogens_min_max_mass(self, h): From 2996d0dd6fd06b9225ebeec04a9bac07c5e67ac2 Mon Sep 17 00:00:00 2001 From: lunamorrow Date: Thu, 25 Apr 2024 15:15:51 +1000 Subject: [PATCH 10/10] revised changes, still needs some adjustments and docstring example updates --- package/CHANGELOG | 2 +- .../analysis/hydrogenbonds/hbond_analysis.py | 76 +++++-------------- 2 files changed, 21 insertions(+), 57 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 44a246745d5..3633a6712c3 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -14,7 +14,7 @@ The rules for this file: ------------------------------------------------------------------------------- -03/27/23 lunamorrow +04/25/24 lunamorrow * 2.8.1 diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index 80d273a61d4..4118132ad9d 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -108,6 +108,10 @@ Alternatively, :attr:`hydrogens_sel` and :attr:`acceptors_sel` may be generated via the :attr:`guess_hydrogens` and :attr:`guess_acceptors`. This selection strings may then be modified prior to calling :attr:`run`, or a subset of +the universe may be used to guess the atoms. For example, find hydrogens and acceptors belonging to a protein:: + +Alternatively, :attr:`hydrogens_ag` and :attr:`acceptors_ag` may be generated via the :attr:`guess_hydrogens` and +:attr:`guess_acceptors`. This selection strings may then be modified prior to calling :attr:`run`, or a subset of the universe may be used to guess the atoms. For example, find hydrogens and acceptors belonging to a protein:: import MDAnalysis @@ -406,9 +410,9 @@ def guess_hydrogens(self, Returns ------- - potential_hydrogens: str - String containing the :attr:`resname` and :attr:`name` of all - hydrogen atoms potentially capable of forming hydrogen bonds. + potential_hydrogens: AtomGroup + AtomGroup containing all of the hydrogen atoms potentially + capable of forming hydrogen bonds. Notes ----- @@ -424,14 +428,11 @@ def guess_hydrogens(self, If :attr:`hydrogens_sel` is `None`, this function is called to guess the selection. - Alternatively, this function may be used to quickly generate a - :class:`str` of potential hydrogen atoms involved in hydrogen bonding. - This str may then be modified before being used to set the attribute - :attr:`hydrogens_sel`. - .. versionchanged:: 2.4.0 Added ability to use atom types + .. versionchanged:: 3.0.0 + Returns AtomGroup instead of selection string """ @@ -439,13 +440,6 @@ def guess_hydrogens(self, raise ValueError("min_mass is higher than (or equal to) max_mass") ag = self.u.select_atoms(select) - # hydrogens_ag = ag[ - # np.logical_and.reduce(( - # ag.masses < max_mass, - # ag.charges > min_charge, - # ag.masses > min_mass, - # )) - # ] hydrogens_ag = ag.select_atoms(f"prop charge > {min_charge} and prop \ mass > {min_mass} and prop mass < \ {max_mass}", @@ -469,9 +463,9 @@ def guess_donors(self, select='all', max_charge=-0.5): Returns ------- - potential_donors: str - String containing the :attr:`resname` and :attr:`name` of all atoms - that are potentially capable of forming hydrogen bonds. + potential_donors: AtomGroup + AtomGroup containing all of the atoms potentially capable of + forming hydrogen bonds. Notes ----- @@ -488,14 +482,11 @@ def guess_donors(self, select='all', max_charge=-0.5): have bonding information, this function is called to guess the selection. - Alternatively, this function may be used to quickly generate a - :class:`str` of potential donor atoms involved in hydrogen bonding. - This :class:`str` may then be modified before being used to set the - attribute :attr:`donors_sel`. - .. versionchanged:: 2.4.0 Added ability to use atom types + .. versionchanged:: 3.0.0 + Returns AtomGroup instead of selection string """ @@ -526,8 +517,6 @@ def guess_donors(self, select='all', max_charge=-0.5): ) ) - #donors_ag = donors_ag[donors_ag.charges < max_charge] - donors_ag = donors_ag.select_atoms(f"prop charge < {max_charge}", updating=self.update_selections) @@ -550,9 +539,9 @@ def guess_acceptors(self, select='all', max_charge=-0.5): Returns ------- - potential_acceptors: str - String containing the :attr:`resname` and :attr:`name` of all atoms - that potentially capable of forming hydrogen bonds. + potential_acceptors: AtomGroup + AtomGroup containing all of the atoms potentially capable of + forming hydrogen bonds. Notes ----- @@ -567,14 +556,11 @@ def guess_acceptors(self, select='all', max_charge=-0.5): If :attr:`acceptors_sel` is `None`, this function is called to guess the selection. - Alternatively, this function may be used to quickly generate a - :class:`str` of potential acceptor atoms involved in hydrogen bonding. - This :class:`str` may then be modified before being used to set the - attribute :attr:`acceptors_sel`. - .. versionchanged:: 2.4.0 Added ability to use atom types + .. versionchanged:: 3.0.0 + Returns AtomGroup instead of selection string """ @@ -658,8 +644,6 @@ def _get_dh_pairs(self): hydrogens = self.u.select_atoms(self.hydrogens_sel, updating=self.update_selections) donors = self.u.select_atoms(self.donors_sel) - # hydrogens = self.guess_hydrogens() - # donors = self.guess_donors() donors_indices, hydrogen_indices = capped_distance( donors.positions, hydrogens.positions, @@ -717,14 +701,8 @@ def _filter_atoms(self, donors, acceptors): def _prepare(self): self.results.hbonds = [[], [], [], [], [], []] - # Set atom selections if they have not been provided - # if self.acceptors_sel is None: - # self.acceptors_sel = self.guess_acceptors() - # if self.hydrogens_sel is None: - # self.hydrogens_sel = self.guess_hydrogens() - # Select atom groups - if self.acceptors_sel == None: + if self.acceptors_sel is None: self._acceptors = self.guess_acceptors() else: self._acceptors = self.u.select_atoms(self.acceptors_sel, @@ -735,14 +713,6 @@ def _single_frame(self): box = self._ts.dimensions - # Update donor-hydrogen pairs if necessary - if self.update_selections: - self._donors, self._hydrogens = self._get_dh_pairs() - - # find D and A within cutoff distance of one another - # min_cutoff = 1.0 as an atom cannot form a hydrogen bond with itself - print(f"D POS: {self._donors.positions}") - print(f"A POS: {self._acceptors.positions}") d_a_indices, d_a_distances = capped_distance( self._donors.positions, self._acceptors.positions, @@ -761,9 +731,6 @@ def _single_frame(self): # Remove D-A pairs more than d_a_cutoff away from one another tmp_donors = self._donors[d_a_indices.T[0]] - print(f"TMP DONORS INITIAL: {tmp_donors}") - print(f"DONORS: {self._donors}") - print(f"DA INDICIES: {d_a_indices.T[0]}") tmp_hydrogens = self._hydrogens[d_a_indices.T[0]] tmp_acceptors = self._acceptors[d_a_indices.T[1]] @@ -796,9 +763,6 @@ def _single_frame(self): # Retrieve atoms, distances and angles of hydrogen bonds hbond_donors = tmp_donors[hbond_indices] - print(f"HBOND DONORS: {hbond_donors}") - print(f"TMP DONORS: {tmp_donors}") - print(f"HBOND INDICIES: {hbond_indices}") hbond_hydrogens = tmp_hydrogens[hbond_indices] hbond_acceptors = tmp_acceptors[hbond_indices] hbond_distances = d_a_distances[hbond_indices]