diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index d4b941614b..2e0d8543d3 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -247,7 +247,9 @@ from MDAnalysis.lib.correlations import autocorrelation, correct_intermittency from MDAnalysis.exceptions import NoDataError from MDAnalysis.core.groups import AtomGroup -from MDAnalysis.analysis.hydrogenbonds.hbond_autocorrel import find_hydrogen_donors +from MDAnalysis.analysis.hydrogenbonds.hbond_autocorrel import ( + find_hydrogen_donors, +) from ...due import due, Doi @@ -255,10 +257,12 @@ logger = logging.getLogger(__name__) -due.cite(Doi("10.1039/C9CP01532A"), - description="Hydrogen bond analysis implementation", - path="MDAnalysis.analysis.hydrogenbonds.hbond_analysis", - cite_module=True) +due.cite( + Doi("10.1039/C9CP01532A"), + description="Hydrogen bond analysis implementation", + path="MDAnalysis.analysis.hydrogenbonds.hbond_analysis", + cite_module=True, +) del Doi @@ -272,18 +276,29 @@ class HydrogenBondAnalysis(AnalysisBase): @classmethod def get_supported_backends(cls): - return ('serial', 'multiprocessing', 'dask',) + return ( + "serial", + "multiprocessing", + "dask", + ) - def __init__(self, universe, - donors_sel=None, hydrogens_sel=None, acceptors_sel=None, - between=None, d_h_cutoff=1.2, - d_a_cutoff=3.0, d_h_a_angle_cutoff=150, - update_selections=True): + def __init__( + self, + universe, + donors_sel=None, + hydrogens_sel=None, + acceptors_sel=None, + between=None, + d_h_cutoff=1.2, + d_a_cutoff=3.0, + d_h_a_angle_cutoff=150, + update_selections=True, + ): """Set up atom selections and geometric criteria for finding hydrogen bonds in a Universe. Hydrogen bond selections with `donors_sel` , `hydrogens_sel`, and - `acceptors_sel` may be achieved with either a *resname*, atom *name* + `acceptors_sel` may be achieved with either a *resname*, atom *name* combination, or when those are absent, with atom *type* selections. Parameters @@ -340,7 +355,7 @@ def __init__(self, universe, .. versionadded:: 2.0.0 Added `between` keyword .. versionchanged:: 2.4.0 - Added use of atom types in selection strings for hydrogen atoms, + Added use of atom types in selection strings for hydrogen atoms, bond donors, or bond acceptors .. versionchanged:: 2.8.0 Introduced :meth:`get_supported_backends` allowing for parallel execution on @@ -354,15 +369,26 @@ def __init__(self, universe, self.u = universe self._trajectory = self.u.trajectory - self._donors_sel = donors_sel.strip() if donors_sel is not None else donors_sel - self._hydrogens_sel = hydrogens_sel.strip() if hydrogens_sel is not None else hydrogens_sel - self._acceptors_sel = acceptors_sel.strip() if acceptors_sel is not None else acceptors_sel + self._donors_sel = ( + donors_sel.strip() if donors_sel is not None else donors_sel + ) + self._hydrogens_sel = ( + hydrogens_sel.strip() + if hydrogens_sel is not None + else hydrogens_sel + ) + self._acceptors_sel = ( + acceptors_sel.strip() + if acceptors_sel is not None + else acceptors_sel + ) - msg = ("{} is an empty selection string - no hydrogen bonds will " - "be found. This may be intended, but please check your " - "selection." - ) - for sel in ['donors_sel', 'hydrogens_sel', 'acceptors_sel']: + msg = ( + "{} is an empty selection string - no hydrogen bonds will " + "be found. This may be intended, but please check your " + "selection." + ) + for sel in ["donors_sel", "hydrogens_sel", "acceptors_sel"]: val = getattr(self, sel) if isinstance(val, str) and not val: warnings.warn(msg.format(sel)) @@ -380,7 +406,7 @@ def __init__(self, universe, between_ags.append( [ self.u.select_atoms(group1, updating=False), - self.u.select_atoms(group2, updating=False) + self.u.select_atoms(group2, updating=False), ] ) @@ -388,7 +414,6 @@ def __init__(self, universe, else: self.between_ags = None - self.d_h_cutoff = d_h_cutoff self.d_a_cutoff = d_a_cutoff self.d_h_a_angle = d_h_a_angle_cutoff @@ -403,23 +428,21 @@ def __init__(self, universe, self._hydrogens_sel = self.guess_hydrogens() # Select atom groups - self._acceptors = self.u.select_atoms(self.acceptors_sel, - updating=self.update_selections) + self._acceptors = self.u.select_atoms( + self.acceptors_sel, updating=self.update_selections + ) self._donors, self._hydrogens = self._get_dh_pairs() - def guess_hydrogens(self, - select='all', - max_mass=1.1, - min_charge=0.3, - min_mass=0.9 - ): + def guess_hydrogens( + self, select="all", max_mass=1.1, min_charge=0.3, min_mass=0.9 + ): """Guesses which hydrogen atoms should be used in the analysis. Parameters ---------- select: str (optional) - :ref:`Selection string ` for atom group - from which hydrogens will be identified. (e.g., ``(resname X and + :ref:`Selection string ` for atom group + from which hydrogens will be identified. (e.g., ``(resname X and name H1)`` or ``type 2``) max_mass: float (optional) The mass of a hydrogen atom must be less than this value. @@ -431,24 +454,24 @@ def guess_hydrogens(self, Returns ------- potential_hydrogens: str - String containing the :attr:`resname` and :attr:`name` of all + String containing the :attr:`resname` and :attr:`name` of all hydrogen atoms potentially capable of forming hydrogen bonds. Notes ----- - Hydrogen selections may be achieved with either a resname, atom + Hydrogen selections may be achieved with either a resname, atom name combination, or when those are absent, atom types. This function makes use of atomic masses and atomic charges to identify - which atoms are hydrogen atoms that are capable of participating in - hydrogen bonding. If an atom has a mass less than :attr:`max_mass` and - an atomic charge greater than :attr:`min_charge` then it is considered + which atoms are hydrogen atoms that are capable of participating in + hydrogen bonding. If an atom has a mass less than :attr:`max_mass` and + an atomic charge greater than :attr:`min_charge` then it is considered capable of participating in hydrogen bonds. - If :attr:`hydrogens_sel` is `None`, this function is called to guess + If :attr:`hydrogens_sel` is `None`, this function is called to guess the selection. - Alternatively, this function may be used to quickly generate a + 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`. @@ -464,25 +487,27 @@ def guess_hydrogens(self, 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, - )) + np.logical_and.reduce( + ( + ag.masses < max_mass, + ag.charges > min_charge, + ag.masses > min_mass, + ) + ) ] return self._group_categories(hydrogens_ag) - def guess_donors(self, select='all', max_charge=-0.5): + def guess_donors(self, select="all", max_charge=-0.5): """Guesses which atoms could be considered donors in the analysis. Only - use if the universe topology does not contain bonding information, + use if the universe topology does not contain bonding information, otherwise donor-hydrogen pairs may be incorrectly assigned. Parameters ---------- select: str (optional) - :ref:`Selection string ` for atom group - from which donors will be identified. (e.g., ``(resname X and name + :ref:`Selection string ` for atom group + from which donors will be identified. (e.g., ``(resname X and name O1)`` or ``type 2``) max_charge: float (optional) The charge of a donor atom must be less than this value. @@ -490,27 +515,27 @@ 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 + String containing the :attr:`resname` and :attr:`name` of all atoms that are potentially capable of forming hydrogen bonds. Notes ----- - Donor selections may be achieved with either a resname, atom + Donor selections may be achieved with either a resname, atom name combination, or when those are absent, atom types. - This function makes use of and atomic charges to identify which atoms - could be considered donor atoms in the hydrogen bond analysis. If an - atom has an atomic charge less than :attr:`max_charge`, and it is - within :attr:`d_h_cutoff` of a hydrogen atom, then it is considered + This function makes use of and atomic charges to identify which atoms + could be considered donor atoms in the hydrogen bond analysis. If an + atom has an atomic charge less than :attr:`max_charge`, and it is + within :attr:`d_h_cutoff` of a hydrogen atom, then it is considered capable of participating in hydrogen bonds. - If :attr:`donors_sel` is `None`, and the universe topology does not - have bonding information, this function is called to guess the + If :attr:`donors_sel` is `None`, and the universe topology does not + 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 + 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`. @@ -520,7 +545,7 @@ 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 + # Use a new variable `hydrogens_sel` so that we do not set # `self.hydrogens_sel` if it is currently `None` if self.hydrogens_sel is None: hydrogens_sel = self.guess_hydrogens() @@ -532,8 +557,10 @@ def guess_donors(self, select='all', max_charge=-0.5): # 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 (hasattr(self.u._topology, 'bonds') - and len(self.u._topology.bonds.values) != 0): + if ( + hasattr(self.u._topology, "bonds") + and len(self.u._topology.bonds.values) != 0 + ): donors_ag = find_hydrogen_donors(hydrogens_ag) donors_ag = donors_ag.intersection(self.u.select_atoms(select)) else: @@ -541,7 +568,7 @@ def guess_donors(self, select='all', max_charge=-0.5): "({donors_sel}) and around {d_h_cutoff} {hydrogens_sel}".format( donors_sel=select, d_h_cutoff=self.d_h_cutoff, - hydrogens_sel=hydrogens_sel + hydrogens_sel=hydrogens_sel, ) ) @@ -549,17 +576,17 @@ def guess_donors(self, select='all', max_charge=-0.5): return self._group_categories(donors_ag) - def guess_acceptors(self, select='all', max_charge=-0.5): + def guess_acceptors(self, select="all", max_charge=-0.5): """Guesses which atoms could be considered acceptors in the analysis. - Acceptor selections may be achieved with either a resname, atom + Acceptor selections may be achieved with either a resname, atom name combination, or when those are absent, atom types. Parameters ---------- select: str (optional) :ref:`Selection string ` for atom group - from which acceptors will be identified. (e.g., ``(resname X and + from which acceptors will be identified. (e.g., ``(resname X and name O1)`` or ``type 2``) max_charge: float (optional) The charge of an acceptor atom must be less than this value. @@ -567,25 +594,25 @@ 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 + String containing the :attr:`resname` and :attr:`name` of all atoms that potentially capable of forming hydrogen bonds. Notes ----- - Acceptor selections may be achieved with either a resname, atom + Acceptor selections may be achieved with either a resname, atom name combination, or when those are absent, atom types. - This function makes use of and atomic charges to identify which atoms - could be considered acceptor atoms in the hydrogen bond analysis. If - an atom has an atomic charge less than :attr:`max_charge` then it is + This function makes use of and atomic charges to identify which atoms + could be considered acceptor atoms in the hydrogen bond analysis. If + an atom has an atomic charge less than :attr:`max_charge` then it is considered capable of participating in hydrogen bonds. - If :attr:`acceptors_sel` is `None`, this function is called to guess + 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 + 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`. @@ -601,14 +628,14 @@ def guess_acceptors(self, select='all', max_charge=-0.5): @staticmethod def _group_categories(group): - """ Find categories according to universe constraints - + """Find categories according to universe constraints + Parameters ---------- group : AtomGroup - AtomGroups corresponding to either hydrogen bond acceptors, + AtomGroups corresponding to either hydrogen bond acceptors, donors, or hydrogen atoms that meet their respective charge - and mass constraints. + and mass constraints. Returns ------- @@ -621,16 +648,16 @@ def _group_categories(group): """ if hasattr(group, "resnames") and hasattr(group, "names"): - group_list = np.unique([ - '(resname {} and name {})'.format(r, - p) for r, p in zip(group.resnames, group.names) - ]) - else: group_list = np.unique( [ - 'type {}'.format(tp) for tp in group.types + "(resname {} and name {})".format(r, p) + for r, p in zip(group.resnames, group.names) ] ) + else: + group_list = np.unique( + ["type {}".format(tp) for tp in group.types] + ) return " or ".join(group_list) @@ -640,7 +667,7 @@ def _get_dh_pairs(self): Returns ------- donors, hydrogens: AtomGroup, AtomGroup - AtomGroups corresponding to all donors and all hydrogens. + AtomGroups corresponding to all donors and all hydrogens. AtomGroups are ordered such that, if zipped, will produce a list of donor-hydrogen pairs. """ @@ -651,15 +678,23 @@ def _get_dh_pairs(self): # 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.') + 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 \ + 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: @@ -671,7 +706,7 @@ def _get_dh_pairs(self): hydrogens.positions, max_cutoff=self.d_h_cutoff, box=self.u.dimensions, - return_distances=False + return_distances=False, ).T donors = donors[donors_indices] @@ -682,20 +717,20 @@ def _get_dh_pairs(self): def _filter_atoms(self, donors, acceptors): """Create a mask to filter donor, hydrogen and acceptor atoms. - This can be used to consider only hydrogen bonds between two or more - specified groups. + This can be used to consider only hydrogen bonds between two or more + specified groups. - Groups are specified with the `between` keyword when creating the - HydrogenBondAnalysis object. + Groups are specified with the `between` keyword when creating the + HydrogenBondAnalysis object. - Returns - ------- - mask: np.ndarray + Returns + ------- + mask: np.ndarray - .. versionchanged:: 2.5.0 - Change return value to a mask instead of separate AtomGroups. -`` + .. versionchanged:: 2.5.0 + Change return value to a mask instead of separate AtomGroups. + `` """ mask = np.full(donors.n_atoms, fill_value=False) @@ -703,27 +738,25 @@ def _filter_atoms(self, donors, acceptors): # Find donors in G1 and acceptors in G2 mask[ - np.logical_and( - np.isin(donors.indices, group1.indices), - np.isin(acceptors.indices, group2.indices) - ) + np.logical_and( + np.isin(donors.indices, group1.indices), + np.isin(acceptors.indices, group2.indices), + ) ] = True # Find acceptors in G1 and donors in G2 mask[ np.logical_and( np.isin(acceptors.indices, group1.indices), - np.isin(donors.indices, group2.indices) + np.isin(donors.indices, group2.indices), ) ] = True return mask - def _prepare(self): self.results.hbonds = [[], [], [], [], [], []] - def _single_frame(self): box = self._ts.dimensions @@ -770,7 +803,7 @@ def _single_frame(self): tmp_donors.positions, tmp_hydrogens.positions, tmp_acceptors.positions, - box=box + box=box, ) ) hbond_indices = np.where(d_h_a_angles > self.d_h_a_angle)[0] @@ -790,8 +823,9 @@ def _single_frame(self): hbond_angles = d_h_a_angles[hbond_indices] # Store data on hydrogen bonds found at this frame - self.results.hbonds[0].extend(np.full_like(hbond_donors, - self._ts.frame)) + self.results.hbonds[0].extend( + np.full_like(hbond_donors, self._ts.frame) + ) self.results.hbonds[1].extend(hbond_donors.indices) self.results.hbonds[2].extend(hbond_hydrogens.indices) self.results.hbonds[3].extend(hbond_acceptors.indices) @@ -803,13 +837,15 @@ def _conclude(self): self.results.hbonds = np.asarray(self.results.hbonds).T def _get_aggregator(self): - return ResultsGroup(lookup={'hbonds': ResultsGroup.ndarray_hstack}) + return ResultsGroup(lookup={"hbonds": ResultsGroup.ndarray_hstack}) @property def hbonds(self): - wmsg = ("The `hbonds` attribute was deprecated in MDAnalysis 2.0.0 " - "and will be removed in MDAnalysis 3.0.0. Please use " - "`results.hbonds` instead.") + wmsg = ( + "The `hbonds` attribute was deprecated in MDAnalysis 2.0.0 " + "and will be removed in MDAnalysis 3.0.0. Please use " + "`results.hbonds` instead." + ) warnings.warn(wmsg, DeprecationWarning) return self.results.hbonds @@ -886,17 +922,16 @@ def lifetime(self, tau_max=20, window_step=1, intermittency=0): # [set(superset(x1,x2), superset(x3,x4)), ..] found_hydrogen_bonds = [set() for _ in self.frames] for frame_index, frame in enumerate(self.frames): - for hbond in self.results.hbonds[self.results.hbonds[:, 0] == frame]: + for hbond in self.results.hbonds[ + self.results.hbonds[:, 0] == frame + ]: found_hydrogen_bonds[frame_index].add(frozenset(hbond[2:4])) intermittent_hbonds = correct_intermittency( - found_hydrogen_bonds, - intermittency=intermittency + found_hydrogen_bonds, intermittency=intermittency ) tau_timeseries, timeseries, timeseries_data = autocorrelation( - intermittent_hbonds, - tau_max, - window_step=window_step + intermittent_hbonds, tau_max, window_step=window_step ) return np.vstack([tau_timeseries, timeseries]) @@ -912,7 +947,9 @@ def count_by_time(self): the number of hydrogen bonds over time. """ hbond_frames = self.results.hbonds[:, 0].astype(int) - frame_unique, frame_counts = np.unique(hbond_frames, return_counts=True) + frame_unique, frame_counts = np.unique( + hbond_frames, return_counts=True + ) frame_min, frame_max = self.frames.min(), self.frames.max() counts = np.zeros(frame_max - frame_min + 1, dtype=int) @@ -925,13 +962,13 @@ def count_by_type(self): Returns ------- counts : numpy.ndarray - Each row of the array contains the donor resname, donor atom type, - acceptor resname, acceptor atom type and the total number of times + Each row of the array contains the donor resname, donor atom type, + acceptor resname, acceptor atom type and the total number of times the hydrogen bond was found. Note ---- - Unique hydrogen bonds are determined through a consideration of the + Unique hydrogen bonds are determined through a consideration of the resname and atom type of the donor and acceptor atoms in a hydrogen bond. """ @@ -947,11 +984,13 @@ def count_by_type(self): tmp_hbonds = np.array([d_res, d.types, a_res, a.types], dtype=str).T hbond_type, type_counts = np.unique( - tmp_hbonds, axis=0, return_counts=True) + tmp_hbonds, axis=0, return_counts=True + ) hbond_type_list = [] for hb_type, hb_count in zip(hbond_type, type_counts): - hbond_type_list.append([":".join(hb_type[:2]), - ":".join(hb_type[2:4]), hb_count]) + hbond_type_list.append( + [":".join(hb_type[:2]), ":".join(hb_type[2:4]), hb_count] + ) return np.array(hbond_type_list) @@ -975,12 +1014,14 @@ def count_by_ids(self): a = self.u.atoms[self.results.hbonds[:, 3].astype(np.intp)] tmp_hbonds = np.array([d.ids, h.ids, a.ids]).T - hbond_ids, ids_counts = np.unique(tmp_hbonds, axis=0, - return_counts=True) + hbond_ids, ids_counts = np.unique( + tmp_hbonds, axis=0, return_counts=True + ) # Find unique hbonds and sort rows so that most frequent observed bonds are at the top of the array - unique_hbonds = np.concatenate((hbond_ids, ids_counts[:, None]), - axis=1) + unique_hbonds = np.concatenate( + (hbond_ids, ids_counts[:, None]), axis=1 + ) unique_hbonds = unique_hbonds[unique_hbonds[:, 3].argsort()[::-1]] return unique_hbonds @@ -988,7 +1029,7 @@ def count_by_ids(self): @property def donors_sel(self): """Selection string for the hydrogen bond donor atoms. - + .. versionadded:: 2.10.0 """ return self._donors_sel @@ -1001,7 +1042,7 @@ def donors_sel(self, value): @property def hydrogens_sel(self): """Selection string for the hydrogen bond hydrogen atoms. - + .. versionadded:: 2.10.0 """ return self._hydrogens_sel @@ -1016,7 +1057,7 @@ def hydrogens_sel(self, value): @property def acceptors_sel(self): """Selection string for the hydrogen bond acceptor atoms. - + .. versionadded:: 2.10.0 """ return self._acceptors_sel @@ -1026,5 +1067,6 @@ def acceptors_sel(self, value): self._acceptors_sel = value if self._acceptors_sel is None: self._acceptors_sel = self.guess_acceptors() - self._acceptors = self.u.select_atoms(self._acceptors_sel, - updating=self.update_selections) \ No newline at end of file + self._acceptors = self.u.select_atoms( + self._acceptors_sel, updating=self.update_selections + ) diff --git a/package/MDAnalysis/coordinates/DCD.py b/package/MDAnalysis/coordinates/DCD.py index 672d5323b5..3eecc6fa13 100644 --- a/package/MDAnalysis/coordinates/DCD.py +++ b/package/MDAnalysis/coordinates/DCD.py @@ -63,7 +63,9 @@ import types import warnings -from .. import units as mdaunits # use mdaunits instead of units to avoid a clash +from .. import ( + units as mdaunits, +) # use mdaunits instead of units to avoid a clash from . import base, core from ..lib.formats.libdcd import DCDFile from ..lib.mdamath import triclinic_box @@ -123,9 +125,10 @@ class DCDReader(base.ReaderBase): .. _wiki: https://github.com/MDAnalysis/mdanalysis/wiki/FileFormats#dcd """ - format = 'DCD' - flavor = 'CHARMM' - units = {'time': 'AKMA', 'length': 'Angstrom'} + + format = "DCD" + flavor = "CHARMM" + units = {"time": "AKMA", "length": "Angstrom"} @store_init_arguments def __init__(self, filename, convert_units=True, dt=None, **kwargs): @@ -146,17 +149,19 @@ def __init__(self, filename, convert_units=True, dt=None, **kwargs): Changed to use libdcd.pyx library and removed the correl function """ super(DCDReader, self).__init__( - filename, convert_units=convert_units, **kwargs) + filename, convert_units=convert_units, **kwargs + ) self._file = DCDFile(self.filename) - self.n_atoms = self._file.header['natoms'] + self.n_atoms = self._file.header["natoms"] - delta = mdaunits.convert(self._file.header['delta'], - self.units['time'], 'ps') + delta = mdaunits.convert( + self._file.header["delta"], self.units["time"], "ps" + ) if dt is None: - dt = delta * self._file.header['nsavc'] - self.skip_timestep = self._file.header['nsavc'] + dt = delta * self._file.header["nsavc"] + self.skip_timestep = self._file.header["nsavc"] - self._ts_kwargs['dt'] = dt + self._ts_kwargs["dt"] = dt self.ts = self._Timestep(self.n_atoms, **self._ts_kwargs) frame = self._file.read() # reset trajectory @@ -168,18 +173,20 @@ def __init__(self, filename, convert_units=True, dt=None, **kwargs): self.ts = self._frame_to_ts(frame, self.ts) # these should only be initialized once self.ts.dt = dt - warnings.warn("DCDReader currently makes independent timesteps" - " by copying self.ts while other readers update" - " self.ts inplace. This behavior will be changed in" - " 3.0 to be the same as other readers. Read more at" - " https://github.com/MDAnalysis/mdanalysis/issues/3889" - " to learn if this change in behavior might affect you.", - category=DeprecationWarning) + warnings.warn( + "DCDReader currently makes independent timesteps" + " by copying self.ts while other readers update" + " self.ts inplace. This behavior will be changed in" + " 3.0 to be the same as other readers. Read more at" + " https://github.com/MDAnalysis/mdanalysis/issues/3889" + " to learn if this change in behavior might affect you.", + category=DeprecationWarning, + ) @staticmethod def parse_n_atoms(filename, **kwargs): with DCDFile(filename) as f: - n_atoms = f.header['natoms'] + n_atoms = f.header["natoms"] return n_atoms def close(self): @@ -196,7 +203,7 @@ def _reopen(self): self.ts.frame = 0 self._frame = -1 self._file.close() - self._file.open('r') + self._file.open("r") def _read_frame(self, i): """read frame i""" @@ -207,9 +214,9 @@ def _read_frame(self, i): def _read_next_timestep(self, ts=None): """copy next frame into timestep""" if self._frame == self.n_frames - 1: - raise IOError('trying to go over trajectory limit') + raise IOError("trying to go over trajectory limit") if ts is None: - #TODO remove copying the ts in 3.0 + # TODO remove copying the ts in 3.0 ts = self.ts.copy() frame = self._file.read() self._frame += 1 @@ -226,21 +233,27 @@ def Writer(self, filename, n_atoms=None, **kwargs): n_atoms=n_atoms, dt=self.ts.dt, convert_units=self.convert_units, - **kwargs) + **kwargs, + ) def _frame_to_ts(self, frame, ts): """convert a dcd-frame to a :class:`TimeStep`""" ts.frame = self._frame - ts.time = (ts.frame + self._file.header['istart']/self._file.header['nsavc']) * self.ts.dt - ts.data['step'] = self._file.tell() + ts.time = ( + ts.frame + self._file.header["istart"] / self._file.header["nsavc"] + ) * self.ts.dt + ts.data["step"] = self._file.tell() # The original unitcell is read as ``[A, gamma, B, beta, alpha, C]`` _ts_order = [0, 2, 5, 4, 3, 1] uc = np.take(frame.unitcell, _ts_order) pi_2 = np.pi / 2 - if (-1.0 <= uc[3] <= 1.0) and (-1.0 <= uc[4] <= 1.0) and ( - -1.0 <= uc[5] <= 1.0): + if ( + (-1.0 <= uc[3] <= 1.0) + and (-1.0 <= uc[4] <= 1.0) + and (-1.0 <= uc[5] <= 1.0) + ): # This file was generated by Charmm, or by NAMD > 2.5, with the # angle cosines of the periodic cell angles written to the DCD # file. This formulation improves rounding behavior for orthogonal @@ -251,7 +264,7 @@ def _frame_to_ts(self, frame, ts): uc[4] = 90.0 - np.arcsin(uc[4]) * 90.0 / pi_2 uc[5] = 90.0 - np.arcsin(uc[5]) * 90.0 / pi_2 # heuristic sanity check: uc = A,B,C,alpha,beta,gamma - elif np.any(uc < 0.) or np.any(uc[3:] > 180.): + elif np.any(uc < 0.0) or np.any(uc[3:] > 180.0): # might be new CHARMM: box matrix vectors H = frame.unitcell.copy() e1, e2, e3 = H[[0, 1, 3]], H[[1, 2, 4]], H[[3, 4, 5]] @@ -273,8 +286,7 @@ def _frame_to_ts(self, frame, ts): @property def dimensions(self): - """unitcell dimensions (*A*, *B*, *C*, *alpha*, *beta*, *gamma*) - """ + """unitcell dimensions (*A*, *B*, *C*, *alpha*, *beta*, *gamma*)""" return self.ts.dimensions @property @@ -282,13 +294,15 @@ def dt(self): """timestep between frames""" return self.ts.dt - def timeseries(self, - asel=None, - atomgroup=None, - start=None, - stop=None, - step=None, - order='afc'): + def timeseries( + self, + asel=None, + atomgroup=None, + start=None, + stop=None, + step=None, + order="afc", + ): """Return a subset of coordinate data for an AtomGroup Parameters @@ -297,7 +311,7 @@ def timeseries(self, The :class:`~MDAnalysis.core.groups.AtomGroup` to read the coordinates from. Defaults to None, in which case the full set of coordinate data is returned. - + .. deprecated:: 2.7.0 asel argument will be renamed to atomgroup in 3.0.0 @@ -332,9 +346,12 @@ def timeseries(self, warnings.warn( "asel argument to timeseries will be renamed to" "'atomgroup' in 3.0, see #3911", - category=DeprecationWarning) + category=DeprecationWarning, + ) if atomgroup: - raise ValueError("Cannot provide both asel and atomgroup kwargs") + raise ValueError( + "Cannot provide both asel and atomgroup kwargs" + ) atomgroup = asel start, stop, step = self.check_slice_indices(start, stop, step) @@ -342,13 +359,15 @@ def timeseries(self, if atomgroup is not None: if len(atomgroup) == 0: raise ValueError( - "Timeseries requires at least one atom to analyze") + "Timeseries requires at least one atom to analyze" + ) atom_numbers = list(atomgroup.indices) else: atom_numbers = list(range(self.n_atoms)) frames = self._file.readframes( - start, stop, step, order=order, indices=atom_numbers) + start, stop, step, order=order, indices=atom_numbers + ) return frames.xyz @@ -360,7 +379,7 @@ class DCDWriter(base.WriterBase): and writes positions in Å and time in AKMA time units. .. warning:: - + Multiple conventions exist for how unit cells are written to DCD files. Until 2.10.0 MDAnalysis followed the old X-PLOR/CHARMM/NAMD (≤2.5) convention of storing ``[A, gamma, B, beta, alpha, C]`` while @@ -381,7 +400,7 @@ class DCDWriter(base.WriterBase): ``[0, 0, 0, 0, 0, 0]``). As this behaviour is poorly defined, it may not match the expectations of other software. - .. versionchanged:: 2.10.0 + .. versionchanged:: 2.10.0 Up to 2.10.0 the :class:`DCDWriter` wrote a unit cell following old NAMD (≤2.5) convention, even though the docs stated that the new NAMD convention was being used. Now the modern NAMD > 2.5 format is @@ -391,21 +410,24 @@ class DCDWriter(base.WriterBase): .. _`Issue #5069`: https://github.com/MDAnalysis/mdanalysis/issues/5069 """ - format = 'DCD' + + format = "DCD" multiframe = True - flavor = 'NAMD' - units = {'time': 'AKMA', 'length': 'Angstrom'} - - def __init__(self, - filename, - n_atoms, - convert_units=True, - step=1, - dt=1, - remarks='', - nsavc=1, - istart=0, - **kwargs): + flavor = "NAMD" + units = {"time": "AKMA", "length": "Angstrom"} + + def __init__( + self, + filename, + n_atoms, + convert_units=True, + step=1, + dt=1, + remarks="", + nsavc=1, + istart=0, + **kwargs, + ): """Parameters ---------- filename : str @@ -442,10 +464,10 @@ def __init__(self, if n_atoms is None: raise ValueError("n_atoms argument is required") self.n_atoms = n_atoms - self._file = DCDFile(self.filename, 'w') + self._file = DCDFile(self.filename, "w") self.step = step self.dt = dt - dt = mdaunits.convert(dt, 'ps', self.units['time']) + dt = mdaunits.convert(dt, "ps", self.units["time"]) delta = float(dt) / nsavc istart = istart if istart is not None else nsavc self._file.write_header( @@ -454,7 +476,8 @@ def __init__(self, nsavc=nsavc, delta=delta, is_periodic=1, - istart=istart) + istart=istart, + ) def _write_next_frame(self, ag): """Write information associated with ``obj`` at current frame into trajectory @@ -488,8 +511,10 @@ def _write_next_frame(self, ag): try: dimensions = ts.dimensions.copy() except AttributeError: - wmsg = ('No dimensions set for current frame, zeroed unitcell ' - 'will be written') + wmsg = ( + "No dimensions set for current frame, zeroed unitcell " + "will be written" + ) warnings.warn(wmsg) dimensions = np.zeros(6) @@ -501,7 +526,7 @@ def _write_next_frame(self, ag): # The DCD unitcell is written as [A, cos(gamma), B, cos(beta), cos(alpha), C] _ts_order = [0, 5, 1, 4, 3, 2] box = np.take(dimensions, _ts_order) - + # Convert angles (indices 1, 3, 4) from degrees to the special cosine format # used by NAMD/VMD: cos(angle) = sin(90 - angle); see Issue 5069 box[[1, 3, 4]] = np.sin(np.deg2rad(90.0 - box[[1, 3, 4]])) diff --git a/package/MDAnalysis/coordinates/DLPoly.py b/package/MDAnalysis/coordinates/DLPoly.py index 71297f978f..503601cc50 100644 --- a/package/MDAnalysis/coordinates/DLPoly.py +++ b/package/MDAnalysis/coordinates/DLPoly.py @@ -35,7 +35,7 @@ from ..lib import util from ..lib.util import cached, store_init_arguments -_DLPOLY_UNITS = {'length': 'Angstrom', 'velocity': 'Angstrom/ps', 'time': 'ps'} +_DLPOLY_UNITS = {"length": "Angstrom", "velocity": "Angstrom/ps", "time": "ps"} class ConfigReader(base.SingleFrameReaderBase): @@ -47,13 +47,14 @@ class ConfigReader(base.SingleFrameReaderBase): coordinates, velocities, and forces are no longer stored in 'F' memory layout, instead now using the numpy default of 'C'. """ - format = 'CONFIG' + + format = "CONFIG" units = _DLPOLY_UNITS def _read_first_frame(self): unitcell = np.zeros((3, 3), dtype=np.float32) - with open(self.filename, 'r') as inf: + with open(self.filename, "r") as inf: self.title = inf.readline().strip() levcfg, imcon, megatm = np.int64(inf.readline().split()[:3]) if not imcon == 0: @@ -113,10 +114,12 @@ def _read_first_frame(self): if has_forces: forces = forces[order] - ts = self.ts = self._Timestep(self.n_atoms, - velocities=has_vels, - forces=has_forces, - **self._ts_kwargs) + ts = self.ts = self._Timestep( + self.n_atoms, + velocities=has_vels, + forces=has_forces, + **self._ts_kwargs, + ) ts._pos = coords if has_vels: ts._velocities = velocities @@ -133,7 +136,8 @@ class HistoryReader(base.ReaderBase): .. versionadded:: 0.11.0 """ - format = 'HISTORY' + + format = "HISTORY" units = _DLPOLY_UNITS @store_init_arguments @@ -142,7 +146,7 @@ def __init__(self, filename, **kwargs): self._cache = {} # "private" file handle - self._file = util.anyopen(self.filename, 'r') + self._file = util.anyopen(self.filename, "r") self.title = self._file.readline().strip() header = np.int64(self._file.readline().split()[:3]) self._levcfg, self._imcon, self.n_atoms = header @@ -157,10 +161,12 @@ def __init__(self, filename, **kwargs): self._has_cell = False self._file.seek(rwnd) - self.ts = self._Timestep(self.n_atoms, - velocities=self._has_vels, - forces=self._has_forces, - **self._ts_kwargs) + self.ts = self._Timestep( + self.n_atoms, + velocities=self._has_vels, + forces=self._has_forces, + **self._ts_kwargs, + ) self._read_next_timestep() def _read_next_timestep(self, ts=None): @@ -168,7 +174,7 @@ def _read_next_timestep(self, ts=None): ts = self.ts line = self._file.readline() # timestep line - if not line.startswith('timestep'): + if not line.startswith("timestep"): raise IOError if self._has_cell: @@ -176,7 +182,7 @@ def _read_next_timestep(self, ts=None): unitcell[0] = self._file.readline().split() unitcell[1] = self._file.readline().split() unitcell[2] = self._file.readline().split() - ts.dimensions = core.triclinic_box(*unitcell) + ts.dimensions = core.triclinic_box(*unitcell) # If ids are given, put them in here # and later sort by them @@ -220,17 +226,17 @@ def _read_frame(self, frame): return self._read_next_timestep() @property - @cached('n_frames') + @cached("n_frames") def n_frames(self): # Second line is traj_key, imcom, n_atoms, n_frames, n_records offsets = [] - with open(self.filename, 'r') as f: + with open(self.filename, "r") as f: f.readline() f.readline() position = f.tell() line = f.readline() - while line.startswith('timestep'): + while line.startswith("timestep"): offsets.append(position) if self._has_cell: f.readline() @@ -251,7 +257,7 @@ def n_frames(self): def _reopen(self): self.close() - self._file = open(self.filename, 'r') + self._file = open(self.filename, "r") self._file.readline() # header is 2 lines self._file.readline() self.ts.frame = -1 diff --git a/package/MDAnalysis/coordinates/TRZ.py b/package/MDAnalysis/coordinates/TRZ.py index 82ada24e6a..4e6a7f20e2 100644 --- a/package/MDAnalysis/coordinates/TRZ.py +++ b/package/MDAnalysis/coordinates/TRZ.py @@ -93,7 +93,7 @@ class TRZReader(base.ReaderBase): format = "TRZ" - units = {'time': 'ps', 'length': 'nm', 'velocity': 'nm/ps'} + units = {"time": "ps", "length": "nm", "velocity": "nm/ps"} @store_init_arguments def __init__(self, trzfilename, n_atoms=None, **kwargs): @@ -114,84 +114,95 @@ def __init__(self, trzfilename, n_atoms=None, **kwargs): If `n_atoms` or the number of atoms in the topology file do not match the number of atoms in the trajectory. """ - wmsg = ("The TRZ reader is deprecated and will be removed in " - "MDAnalysis version 3.0.0") + wmsg = ( + "The TRZ reader is deprecated and will be removed in " + "MDAnalysis version 3.0.0" + ) warnings.warn(wmsg, DeprecationWarning) - super(TRZReader, self).__init__(trzfilename, **kwargs) + super(TRZReader, self).__init__(trzfilename, **kwargs) if n_atoms is None: - raise ValueError('TRZReader requires the n_atoms keyword') + raise ValueError("TRZReader requires the n_atoms keyword") - self.trzfile = util.anyopen(self.filename, 'rb') + self.trzfile = util.anyopen(self.filename, "rb") self._cache = dict() self._n_atoms = n_atoms self._read_trz_header() - self.ts = Timestep(self.n_atoms, - velocities=True, - forces=self.has_force, - reader=self, - **self._ts_kwargs) + self.ts = Timestep( + self.n_atoms, + velocities=True, + forces=self.has_force, + reader=self, + **self._ts_kwargs, + ) # structured dtype of a single trajectory frame - readarg = str(n_atoms) + ' angstroms) @@ -245,7 +257,7 @@ def n_atoms(self): return self._n_atoms @property - @cached('n_frames') + @cached("n_frames") def n_frames(self): """Total number of frames in a trajectory""" try: @@ -261,10 +273,17 @@ def _read_trz_n_frames(self, trzfile): """ fsize = os.fstat(trzfile.fileno()).st_size # size of file in bytes - if not (fsize - self._headerdtype.itemsize) % self._dtype.itemsize == 0: - raise IOError("Trajectory has incomplete frames") # check that division is sane + if ( + not (fsize - self._headerdtype.itemsize) % self._dtype.itemsize + == 0 + ): + raise IOError( + "Trajectory has incomplete frames" + ) # check that division is sane - nframes = int((fsize - self._headerdtype.itemsize) / self._dtype.itemsize) # returns long int otherwise + nframes = int( + (fsize - self._headerdtype.itemsize) / self._dtype.itemsize + ) # returns long int otherwise return nframes @@ -294,20 +313,20 @@ def _get_dt(self): self._read_frame(curr_frame) @property - @cached('delta') + @cached("delta") def delta(self): """MD integration timestep""" return self.dt / self.skip_timestep @property - @cached('skip_timestep') + @cached("skip_timestep") def skip_timestep(self): """Timesteps between trajectory frames""" curr_frame = self.ts.frame try: - t0 = self.ts.data['frame'] + t0 = self.ts.data["frame"] self.next() - t1 = self.ts.data['frame'] + t1 = self.ts.data["frame"] skip_timestep = t1 - t0 except StopIteration: return 0 @@ -346,8 +365,12 @@ def _seek(self, nframes): framesize = long(framesize) seeksize = framesize * nframes - nreps = int(seeksize / maxi_l) # number of max seeks we'll have to do - rem = int(seeksize % maxi_l) # amount leftover to do once max seeks done + nreps = int( + seeksize / maxi_l + ) # number of max seeks we'll have to do + rem = int( + seeksize % maxi_l + ) # amount leftover to do once max seeks done for _ in range(nreps): self.trzfile.seek(maxint, 1) @@ -365,13 +388,15 @@ def _reopen(self): def open_trajectory(self): """Open the trajectory file""" if self.trzfile is not None: - raise IOError(errno.EALREADY, 'TRZ file already opened', self.filename) + raise IOError( + errno.EALREADY, "TRZ file already opened", self.filename + ) if not os.path.exists(self.filename): - raise IOError(errno.ENOENT, 'TRZ file not found', self.filename) + raise IOError(errno.ENOENT, "TRZ file not found", self.filename) - self.trzfile = util.anyopen(self.filename, 'rb') + self.trzfile = util.anyopen(self.filename, "rb") - #Reset ts + # Reset ts ts = self.ts ts.frame = -1 @@ -403,12 +428,12 @@ class TRZWriter(base.WriterBase): and will be removed in version 3.0.0. """ - format = 'TRZ' + format = "TRZ" multiframe = True - units = {'time': 'ps', 'length': 'nm', 'velocity': 'nm/ps'} + units = {"time": "ps", "length": "nm", "velocity": "nm/ps"} - def __init__(self, filename, n_atoms, title='TRZ', convert_units=True): + def __init__(self, filename, n_atoms, title="TRZ", convert_units=True): """Create a TRZWriter Parameters @@ -423,8 +448,10 @@ def __init__(self, filename, n_atoms, title='TRZ', convert_units=True): convert_units : bool (optional) units are converted to the MDAnalysis base format; [``True``] """ - wmsg = ("The TRZ writer is deprecated and will be removed in " - "MDAnalysis version 3.0.0") + wmsg = ( + "The TRZ writer is deprecated and will be removed in " + "MDAnalysis version 3.0.0" + ) warnings.warn(wmsg, DeprecationWarning) self.filename = filename @@ -435,65 +462,77 @@ def __init__(self, filename, n_atoms, title='TRZ', convert_units=True): self.n_atoms = n_atoms if len(title) > 80: - raise ValueError("TRZWriter: 'title' must be 80 characters of shorter") + raise ValueError( + "TRZWriter: 'title' must be 80 characters of shorter" + ) self.convert_units = convert_units - self.trzfile = util.anyopen(self.filename, 'wb') + self.trzfile = util.anyopen(self.filename, "wb") self._writeheader(title) - floatsize = str(n_atoms) + 'q') + return self._unpack_value(8, ">q") def unpack_uint64(self): - return self._unpack_value(8, '>Q') + return self._unpack_value(8, ">Q") def unpack_ushort(self): return self.unpack_uint() @@ -114,7 +115,7 @@ def unpack_ushort(self): def unpack_uchar(self): # TPX files prior to gromacs 2020 (tpx version < 119) use unsigned ints # (4 bytes) instead of unsigned chars. - return self._unpack_value(4, '>I') + return self._unpack_value(4, ">I") def do_string(self): """ @@ -138,11 +139,12 @@ class TPXUnpacker2020(TPXUnpacker): gromacs 2020, changes le meaning of some types in the file body (the header keep using the previous implementation of the serializer). """ + @classmethod def from_unpacker(cls, unpacker): new_unpacker = cls(unpacker.get_buffer()) new_unpacker._pos = unpacker.get_position() - if hasattr(unpacker, 'unpack_real'): + if hasattr(unpacker, "unpack_real"): if unpacker.unpack_real == unpacker.unpack_float: new_unpacker.unpack_real = new_unpacker.unpack_float elif unpacker.unpack_real == unpacker.unpack_double: @@ -153,7 +155,7 @@ def from_unpacker(cls, unpacker): def unpack_fstring(self, n): if n < 0: - raise ValueError('Size of fstring cannot be negative.') + raise ValueError("Size of fstring cannot be negative.") start_position = self._pos end_position = self._pos = start_position + n if end_position > len(self._buf): @@ -164,12 +166,12 @@ def unpack_fstring(self, n): def unpack_ushort(self): # The InMemorySerializer implements ushort according to the XDR standard # on the contrary to the IO serializer. - return self._unpack_value(2, '>H') + return self._unpack_value(2, ">H") def unpack_uchar(self): # The InMemorySerializer implements uchar according to the XDR standard # on the contrary to the IO serializer. - return self._unpack_value(1, '>B') + return self._unpack_value(1, ">B") def do_string(self): """ @@ -195,7 +197,9 @@ def do_rvec(data): def ndo_rvec(data, n): """mimic of gmx_fio_ndo_rvec in gromacs""" - return [data.unpack_farray(setting.DIM, data.unpack_real) for i in range(n)] + return [ + data.unpack_farray(setting.DIM, data.unpack_real) for i in range(n) + ] def ndo_ivec(data, n): @@ -221,12 +225,11 @@ def define_unpack_real(prec, data): def read_tpxheader(data): - """this function is now compatible with do_tpxheader in tpxio.cpp - """ + """this function is now compatible with do_tpxheader in tpxio.cpp""" # Last compatibility check with gromacs-2016 ver_str = data.do_string() # version string e.g. VERSION 4.0.5 - if not ver_str.startswith(b'VERSION'): - raise ValueError('Input does not look like a TPR file.') + if not ver_str.startswith(b"VERSION"): + raise ValueError("Input does not look like a TPR file.") precision = data.unpack_int() # e.g. 4 define_unpack_real(precision, data) fileVersion = data.unpack_int() # version of tpx file @@ -239,17 +242,23 @@ def read_tpxheader(data): data.unpack_int() # the value is 8, but haven't found the file_tag = data.do_string() - fileGeneration = data.unpack_int() if fileVersion >= 26 else 0 # generation of tpx file, e.g. 17 + fileGeneration = ( + data.unpack_int() if fileVersion >= 26 else 0 + ) # generation of tpx file, e.g. 17 # Versions before 77 don't have the tag, set it to TPX_TAG_RELEASE file_tag # file_tag is used for comparing with tpx_tag. Only tpr files with a # tpx_tag from a lower or the same version of gromacs code can be parsed by # the tpxio.c - file_tag = data.do_string() if fileVersion >= 81 else setting.TPX_TAG_RELEASE + file_tag = ( + data.do_string() if fileVersion >= 81 else setting.TPX_TAG_RELEASE + ) natoms = data.unpack_int() # total number of atoms - ngtc = data.unpack_int() if fileVersion >= 28 else 0 # number of groups for T-coupling + ngtc = ( + data.unpack_int() if fileVersion >= 28 else 0 + ) # number of groups for T-coupling if fileVersion < 62: # not sure what these two are for. @@ -272,9 +281,24 @@ def read_tpxheader(data): if fileVersion >= setting.tpxv_AddSizeField and fileGeneration >= 27: sizeOfTprBody = data.unpack_int64() - th = obj.TpxHeader(ver_str, precision, fileVersion, fileGeneration, - file_tag, natoms, ngtc, fep_state, lamb, - bIr, bTop, bX, bV, bF, bBox, sizeOfTprBody) + th = obj.TpxHeader( + ver_str, + precision, + fileVersion, + fileGeneration, + file_tag, + natoms, + ngtc, + fep_state, + lamb, + bIr, + bTop, + bX, + bV, + bF, + bBox, + sizeOfTprBody, + ) return th @@ -328,9 +352,11 @@ def do_mtop(data, fver, tpr_resid_from_one=False, precision=4): mb = do_molblock(data) # segment is made to correspond to the molblock as in gromacs, the # naming is kind of arbitrary - molblock = mtop.moltypes[mb.molb_type].name.decode('utf-8') + molblock = mtop.moltypes[mb.molb_type].name.decode("utf-8") segid = f"seg_{i}_{molblock}" - chainID = molblock[14:] if molblock[:14] == "Protein_chain_" else molblock + chainID = ( + molblock[14:] if molblock[:14] == "Protein_chain_" else molblock + ) for j in range(mb.molb_nmol): mt = mtop.moltypes[mb.molb_type] # mt: molecule type for atomkind in mt.atomkinds: @@ -372,17 +398,11 @@ def do_mtop(data, fver, tpr_resid_from_one=False, precision=4): resids += 1 resnames = np.array(resnames, dtype=object) - (residx, new_resids, - (new_resnames, - new_moltypes, - new_molnums, - perres_segids - ) - ) = squash_by(resids, - resnames, - moltypes, - molnums, - segids) + ( + residx, + new_resids, + (new_resnames, new_moltypes, new_molnums, perres_segids), + ) = squash_by(resids, resnames, moltypes, molnums, segids) residueids = Resids(new_resids) residuenames = Resnames(new_resnames) residue_moltypes = Moltypes(new_moltypes) @@ -414,10 +434,12 @@ def do_mtop(data, fver, tpr_resid_from_one=False, precision=4): ) top.add_TopologyAttr(Bonds([bond for bond in bonds if bond])) top.add_TopologyAttr(Angles([angle for angle in angles if angle])) - top.add_TopologyAttr(Dihedrals([dihedral for dihedral in dihedrals - if dihedral])) - top.add_TopologyAttr(Impropers([improper for improper in impropers - if improper])) + top.add_TopologyAttr( + Dihedrals([dihedral for dihedral in dihedrals if dihedral]) + ) + top.add_TopologyAttr( + Impropers([improper for improper in impropers if improper]) + ) if any(elements): elements = Elements(np.array(elements, dtype=object)) @@ -467,7 +489,7 @@ def do_mtop(data, fver, tpr_resid_from_one=False, precision=4): if fver > 58: ngrid = data.unpack_int() grid_spacing = data.unpack_int() - n_elements = grid_spacing ** 2 + n_elements = grid_spacing**2 for i in range(ngrid): for j in range(n_elements): ndo_real(data, 4) @@ -547,9 +569,12 @@ def do_iparams(data, functypes, fver): # Not all elif cases in this function has been used and tested for k, i in enumerate(functypes): if i in [ - setting.F_ANGLES, setting.F_G96ANGLES, - setting.F_BONDS, setting.F_G96BONDS, - setting.F_HARMONIC, setting.F_IDIHS + setting.F_ANGLES, + setting.F_G96ANGLES, + setting.F_BONDS, + setting.F_G96BONDS, + setting.F_HARMONIC, + setting.F_IDIHS, ]: do_harm(data) elif i in [setting.F_RESTRANGLES]: @@ -577,8 +602,10 @@ def do_iparams(data, functypes, fver): data.unpack_real() # restraint.up2B data.unpack_real() # restraint.kB elif i in [ - setting.F_TABBONDS, setting.F_TABBONDSNC, - setting.F_TABANGLES, setting.F_TABDIHS, + setting.F_TABBONDS, + setting.F_TABBONDSNC, + setting.F_TABANGLES, + setting.F_TABDIHS, ]: data.unpack_real() # tab.kA data.unpack_int() # tab.table @@ -604,7 +631,7 @@ def do_iparams(data, functypes, fver): data.unpack_real() # u_b.kUBB elif i in [setting.F_QUARTIC_ANGLES]: data.unpack_real() # qangle.theta - ndo_real(data, 5) # qangle.c + ndo_real(data, 5) # qangle.c elif i in [setting.F_BHAM]: data.unpack_real() # bham.a data.unpack_real() # bham.b @@ -664,8 +691,10 @@ def do_iparams(data, functypes, fver): data.unpack_real() # ljcnb.c12 elif i in [ - setting.F_PIDIHS, setting.F_ANGRES, - setting.F_ANGRESZ, setting.F_PDIHS, + setting.F_PIDIHS, + setting.F_ANGRES, + setting.F_ANGRESZ, + setting.F_PDIHS, ]: data.unpack_real() # pdihs_phiA data.unpack_real() # pdihs_cpA @@ -714,8 +743,8 @@ def do_iparams(data, functypes, fver): do_rvec(data) # posres.fcB elif i in [setting.F_FBPOSRES]: - data.unpack_int() # fbposres.geom - do_rvec(data) # fbposres.pos0 + data.unpack_int() # fbposres.geom + do_rvec(data) # fbposres.pos0 data.unpack_real() # fbposres.r data.unpack_real() # fbposres.k @@ -751,7 +780,11 @@ def do_iparams(data, functypes, fver): data.unpack_real() # vsite.a data.unpack_real() # vsite.b - elif i in [setting.F_VSITE3OUT, setting.F_VSITE4FD, setting.F_VSITE4FDN]: + elif i in [ + setting.F_VSITE3OUT, + setting.F_VSITE4FD, + setting.F_VSITE4FDN, + ]: data.unpack_real() # vsite.a data.unpack_real() # vsite.b data.unpack_real() # vsite.c @@ -794,16 +827,18 @@ def do_moltype(data, symtab, fver): #### start: MDAnalysis specific atomkinds = [] for k, a in enumerate(atoms_obj.atoms): - atomkinds.append(obj.AtomKind( - k, - atoms_obj.atomnames[k], - atoms_obj.type[k], - a.resind, - atoms_obj.resnames[a.resind], - a.m, - a.q, - a.atomnumber, - )) + atomkinds.append( + obj.AtomKind( + k, + atoms_obj.atomnames[k], + atoms_obj.type[k], + a.resind, + atoms_obj.resnames[a.resind], + a.m, + a.q, + a.atomnumber, + ) + ) #### end: MDAnalysis specific # key info: about bonds, angles, dih, improp dih. @@ -819,21 +854,45 @@ def do_moltype(data, symtab, fver): # the following if..elif..else statement needs to be updated as new # type of interactions become interested - if ik_obj.name in ['BONDS', 'G96BONDS', 'MORSE', 'CUBICBONDS', - 'CONNBONDS', 'HARMONIC', 'FENEBONDS', - 'RESTRAINTPOT', 'CONSTR', 'CONSTRNC', - 'TABBONDS', 'TABBONDSNC']: + if ik_obj.name in [ + "BONDS", + "G96BONDS", + "MORSE", + "CUBICBONDS", + "CONNBONDS", + "HARMONIC", + "FENEBONDS", + "RESTRAINTPOT", + "CONSTR", + "CONSTRNC", + "TABBONDS", + "TABBONDSNC", + ]: bonds += list(ik_obj.process(ias)) - elif ik_obj.name in ['ANGLES', 'G96ANGLES', 'CROSS_BOND_BOND', - 'CROSS_BOND_ANGLE', 'UREY_BRADLEY', 'QANGLES', - 'RESTRANGLES', 'TABANGLES', 'LINEAR_ANGLES']: + elif ik_obj.name in [ + "ANGLES", + "G96ANGLES", + "CROSS_BOND_BOND", + "CROSS_BOND_ANGLE", + "UREY_BRADLEY", + "QANGLES", + "RESTRANGLES", + "TABANGLES", + "LINEAR_ANGLES", + ]: angles += list(ik_obj.process(ias)) - elif ik_obj.name in ['PDIHS', 'RBDIHS', 'RESTRDIHS', 'CBTDIHS', - 'FOURDIHS', 'TABDIHS']: + elif ik_obj.name in [ + "PDIHS", + "RBDIHS", + "RESTRDIHS", + "CBTDIHS", + "FOURDIHS", + "TABDIHS", + ]: dihs += list(ik_obj.process(ias)) - elif ik_obj.name in ['IDIHS', 'PIDIHS']: + elif ik_obj.name in ["IDIHS", "PIDIHS"]: impr += list(ik_obj.process(ias)) - elif ik_obj.name == 'SETTLE': + elif ik_obj.name == "SETTLE": # SETTLE interactions are optimized triangular constraints for # water molecules. They should be counted as a pair of bonds # between the oxygen and the hydrogens. In older versions of @@ -922,7 +981,9 @@ def do_ilists(data, fver): nr = [] # number of ilist iatoms = [] # atoms involved in a particular interaction type pos = [] - for j in range(setting.F_NRE): # total number of energies (i.e. interaction types) + for j in range( + setting.F_NRE + ): # total number of energies (i.e. interaction types) bClear = False for k in setting.ftupd: if fver < k[0] and j == k[1]: @@ -952,12 +1013,15 @@ def do_molblock(data): molb_nposres_xA = data.unpack_int() if molb_nposres_xA > 0: ndo_rvec(data, molb_nposres_xA) - molb_nposres_xB = data.unpack_int() # The number of posres coords for top B + molb_nposres_xB = ( + data.unpack_int() + ) # The number of posres coords for top B if molb_nposres_xB > 0: ndo_rvec(data, molb_nposres_xB) - return obj.Molblock(molb_type, molb_nmol, molb_natoms_mol, - molb_nposres_xA, molb_nposres_xB) + return obj.Molblock( + molb_type, molb_nmol, molb_natoms_mol, molb_nposres_xA, molb_nposres_xB + ) def do_block(data): @@ -969,7 +1033,9 @@ def do_block(data): def do_blocka(data): block_nr = data.unpack_int() # No. of atoms with excls - block_nra = data.unpack_int() # total times fo appearance of atoms for excls + block_nra = ( + data.unpack_int() + ) # total times fo appearance of atoms for excls ndo_int(data, block_nr + 1) ndo_int(data, block_nra) return block_nr, block_nra diff --git a/package/pyproject.toml b/package/pyproject.toml index 283980b133..fe6e79da74 100644 --- a/package/pyproject.toml +++ b/package/pyproject.toml @@ -130,16 +130,11 @@ line-length = 79 target-version = ['py311', 'py312', 'py313'] extend-exclude = ''' ( -__pycache__ -| MDAnalysis/coordinates/DCD\.py +__pycache__ | MDAnalysis/coordinates/TRJ\.py | MDAnalysis/coordinates/__init__\.py -| MDAnalysis/analysis/hydrogenbonds/hbond_analysis\.py -| MDAnalysis/coordinates/DLPoly\.py -| MDAnalysis/coordinates/TRZ\.py | MDAnalysis/coordinates/__init__\.py | MDAnalysis/coordinates/base\.py -| MDAnalysis/topology/tpr/utils\.py | MDAnalysis/coordinates/MMCIF\.py | MDAnalysis/topology/MMCIFParser\.py | MDAnalysis/topology/PDBParser\.py