From 224ad04b0af452e120b40f833b9809d393bd0ef1 Mon Sep 17 00:00:00 2001
From: Rocco Meli <r.meli@bluemail.ch>
Date: Wed, 2 Oct 2024 14:38:53 +0200
Subject: [PATCH 1/3] black action

---
 .github/workflows/darkerbot.yaml |  62 --------
 .github/workflows/linters.yaml   |  71 ++-------
 maintainer/ci/darker-outcomes.py | 238 -------------------------------
 3 files changed, 10 insertions(+), 361 deletions(-)
 delete mode 100644 .github/workflows/darkerbot.yaml
 delete mode 100644 maintainer/ci/darker-outcomes.py

diff --git a/.github/workflows/darkerbot.yaml b/.github/workflows/darkerbot.yaml
deleted file mode 100644
index 88d0a085159..00000000000
--- a/.github/workflows/darkerbot.yaml
+++ /dev/null
@@ -1,62 +0,0 @@
-name: Darker PR Bot
-
-on:
-  workflow_run:
-    workflows: [linters]
-    types:
-      - completed
-
-
-concurrency:
-  # Probably overly cautious group naming.
-  # Commits to develop will cancel each other, but PRs will only cancel
-  # commits within the same PR
-  group: "${{ github.ref }}-${{ github.head_ref }}-${{ github.workflow }}"
-  cancel-in-progress: true
-
-
-jobs:
-  darker_bot:
-    if: "github.repository == 'MDAnalysis/mdanalysis'"
-    runs-on: ubuntu-latest
-    timeout-minutes: 5
-    steps:
-      - uses: actions/checkout@v4
-
-      - uses: actions/setup-python@v5
-        with:
-          python-version: "3.10"
-
-      - name: setup_dependencies
-        run: |
-          pip install PyGithub
-
-      - name: 'Download artifact'
-        uses: actions/github-script@v7
-        with:
-          script: |
-            let allArtifacts = await github.rest.actions.listWorkflowRunArtifacts({
-               owner: context.repo.owner,
-               repo: context.repo.repo,
-               run_id: context.payload.workflow_run.id,
-            });
-            let matchArtifact = allArtifacts.data.artifacts.filter((artifact) => {
-              return artifact.name == "darkerlint"
-            })[0];
-            let download = await github.rest.actions.downloadArtifact({
-               owner: context.repo.owner,
-               repo: context.repo.repo,
-               artifact_id: matchArtifact.id,
-               archive_format: 'zip',
-            });
-            let fs = require('fs');
-            fs.writeFileSync(`${process.env.GITHUB_WORKSPACE}/darker_results.zip`, Buffer.from(download.data));
-
-      - name: 'Unzip artifact'
-        run: unzip darker_results.zip
-
-      - name: writer-errors
-        env:
-          GITHUB_TOKEN: ${{ github.token }}
-        run: |
-          python maintainer/ci/darker-outcomes.py --json status.json
diff --git a/.github/workflows/linters.yaml b/.github/workflows/linters.yaml
index de5eb68fa86..ebc6225036c 100644
--- a/.github/workflows/linters.yaml
+++ b/.github/workflows/linters.yaml
@@ -16,12 +16,10 @@ defaults:
     shell: bash -l {0}
 
 jobs:
-  darker_lint:
+  black:
     if: "github.repository == 'MDAnalysis/mdanalysis'"
     runs-on: ubuntu-latest
     timeout-minutes: 10
-    permissions:
-      pull-requests: write
     defaults:
       run:
         shell: bash
@@ -35,66 +33,17 @@ jobs:
       with:
         python-version: "3.10"
 
-    - name: darker-main-code
-      id: darker-main-code
-      uses: akaihola/darker@v2.1.1
-      continue-on-error: true
+    - uses: psf/black@stable
       with:
-        version: "~=1.6.1"
-        options: "--check --diff --color"
-        src: "./package/MDAnalysis"
-        revision: "HEAD^"
-        lint: "flake8"
-
-    - name: darker-test-code
-      id: darker-test-code
-      uses: akaihola/darker@v2.1.1
-      continue-on-error: true
+        options: "--check --verbose"
+        src: "./package"
+        version: "~= 24.0"
+    
+    - uses: psf/black@stable
       with:
-        version: "~=1.6.1"
-        options: "--check --diff --color"
-        src: "./testsuite/MDAnalysisTests"
-        revision: "HEAD^"
-        lint: "flake8"
-
-    - name: get-pr-info
-      uses: actions/github-script@v7
-      with:
-        script:
-            const prNumber = context.payload.number;
-            core.exportVariable('PULL_NUMBER', prNumber);
-
-    - name: save-status
-      env:
-        MAIN: ${{ steps.darker-main-code.outcome }}
-        TEST: ${{ steps.darker-test-code.outcome }}
-      shell: python
-      run: |
-        import os
-        import json
-        from pathlib import Path
-
-        Path('./darker_results/').mkdir(exist_ok=True)
-
-        d = {
-            'main_stat': os.environ['MAIN'],
-            'test_stat': os.environ['TEST'],
-            'PR_NUM': os.environ['PULL_NUMBER'],
-            'RUN_ID': os.environ['GITHUB_RUN_ID'],
-        }
-
-        with open('darker_results/status.json', 'w') as f:
-            json.dump(d, f)
-
-    - name: check-json
-      run: cat darker_results/status.json
-
-    - name: upload-status
-      uses: actions/upload-artifact@v4
-      with:
-        name: darkerlint
-        path: darker_results/
-        retention-days: 1
+        options: "--check --verbose"
+        src: "./testsuite"
+        version: "~= 24.0"
 
 
   pylint_check:
diff --git a/maintainer/ci/darker-outcomes.py b/maintainer/ci/darker-outcomes.py
deleted file mode 100644
index 4d6bb7fc5ba..00000000000
--- a/maintainer/ci/darker-outcomes.py
+++ /dev/null
@@ -1,238 +0,0 @@
-#!/usr/bin/env python3
-
-# MIT License
-
-# Copyright (c) 2023 Irfan Alibay
-
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-
-# The above copyright notice and this permission notice shall be included in
-# all copies or substantial portions of the Software.
-
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-# SOFTWARE.
-
-import argparse
-import os
-from urllib import request
-import json
-from github import Github
-
-
-parser = argparse.ArgumentParser(
-    description="Write PR comment for failed darker linting",
-)
-
-
-parser.add_argument(
-    "--json",
-    type=str,
-    help="Input JSON file with status results",
-)
-
-
-def get_pull_request(repo, pr_num):
-    """
-    Simple method to get a PyGithub PR object from a PR number
-
-    Parameters
-    ----------
-    repo : PyGithub.Repository
-      Github Repository API object.
-    pr_num : int
-      Pull request number.
-
-    Returns
-    -------
-    PyGithub.PullRequest
-      Pull Request object corresponding to the input PR number.
-    """
-    # Can get PR directly from PR number
-    return repo.get_pull(pr_num)
-
-
-def get_action_url(repo, pr, run_id, workflow_name, job_name):
-    """
-    Mix PyGithub & Github V3 REST API method to extract the url path to a
-    github actions workflow job corresponding to a given GITHUB_RUN_ID.
-
-    Parameters
-    ----------
-    repo : PyGithub.Repository
-      Github Repository API object.
-    run_id : str
-      Github actions RUN ID as defined by the github actions environment
-      variable `GITHUB_RUN_ID`.
-    workflow_name : str
-      Name of the workflow to extract a job url for.
-    job_name : str
-      Name of the job within the workflow to extract a url for.
-
-
-    Returns
-    -------
-    str
-      URL to github actions workflow job or 'N/A' if a corresponding job name
-      could not be found.
-    """
-    # Accessing get_workflow directly currently fails when passing a name
-    # Instead do the roundabout way by getting list of all workflows
-    linters = [wf for wf in repo.get_workflows()
-               if wf.name == workflow_name][0]
-
-    # Extract the gh action run
-    run = [r for r in linters.get_runs(branch=pr.head.ref)
-           if r.id == int(run_id)][0]
-
-    # The exact job url can't be recovered via the Python API
-    # Switch over to using the REST API instead
-    with request.urlopen(run.jobs_url) as url:
-        data = json.load(url)
-
-    for job in data['jobs']:
-        if job['name'] == job_name:
-            return job['html_url']
-
-    return 'N/A'
-
-
-def bool_outcome(outcome):
-    """
-    Converts github action job status outcome to bool.
-
-    Parameters
-    ----------
-    outcome : str
-      Github action job step outcome message.
-
-    Returns
-    -------
-    bool
-      Whether or not the job step was successful.
-    """
-    return True if (outcome == 'success') else False
-
-
-def gen_message(pr, main_stat, test_stat, action_url):
-    """
-    Generate a user facing message on the status of the darker linting
-    action.
-
-    Parameters
-    ----------
-    pr : PyGithub.PullRequest
-      Pull Request object representing the target for this message.
-    main_stat : str
-      Outcome of darker linting of main package code.
-    test_stat : str
-      Outcome of darker linting of testsuite code.
-    action_url : str
-      URL pointing to darker linting job log.
-
-
-    Returns
-    -------
-    str
-      Message to be posted to PR author.
-    """
-
-    def _format_outcome(stat):
-        if bool_outcome(stat):
-            return "✅ Passed"
-        else:
-            return "⚠️  Possible failure"
-
-    msg = ('### Linter Bot Results:\n\n'
-           f'Hi @{pr.user.login}! Thanks for making this PR. '
-           'We linted your code and found the following: \n\n')
-
-    # If everything is ok
-    if bool_outcome(main_stat) and bool_outcome(test_stat):
-        msg += ('There are currently no issues detected! 🎉')
-    else:
-        msg += ('Some issues were found with the formatting of your code.\n'
-                '| Code Location | Outcome |\n'
-                '| --- | --- |\n'
-                f'| main package | {_format_outcome(main_stat)}|\n'
-                f'| testsuite | {_format_outcome(test_stat)}|\n'
-                '\nPlease have a look at the `darker-main-code` and '
-                '`darker-test-code` steps here for more details: '
-                f'{action_url}\n\n'
-                '---\n'
-                '_**Please note:** The `black` linter is purely '
-                'informational, you can safely ignore these outcomes if '
-                'there are no flake8 failures!_')
-    return msg
-
-
-def post_comment(pr, message, match_string):
-    """
-    Post a comment in a Pull Request.
-
-    If a comment with text matching `match_string` is found in the
-    Pull Request, the comment will be edited.
-
-    Parameters
-    ----------
-    pr : PyGithub.PullRequest
-      Pull Request object representing the target for this message.
-    message : str
-      The message to post as a comment.
-    match_string : str
-      A matching string to recognise if the comment already exists.
-    """
-    # Extract a list of matching comments from PR
-    comments = [comm for comm in pr.get_issue_comments() if match_string in comm.body]
-
-    if len(comments) > 0:
-        # Edit the comment in-place
-        # By default we assume that the bot can write faster than anyone else
-        comments[0].edit(body=message)
-    else:
-        # If the comment doesn't exist, create one
-        pr.create_issue_comment(message)
-
-
-if __name__ == "__main__":
-    args = parser.parse_args()
-
-    git = Github(os.environ['GITHUB_TOKEN'])
-    repo = git.get_repo("MDAnalysis/mdanalysis")
-
-    with open(args.json, 'r') as f:
-        status = json.load(f)
-
-    run_id = status['RUN_ID']
-    print(f"debug run_id: {run_id}")
-
-    # Get Pull Request
-    pr_num = int(status['PR_NUM'])
-    print(f"debug pr_num: {pr_num}")
-    pr = get_pull_request(repo, pr_num)
-
-    # Get the url to the github action job being pointed to
-    action_url = get_action_url(repo, pr, run_id,
-                                workflow_name='linters',
-                                job_name='darker_lint')
-
-    # Get the message you want to post to users
-    with open(args.json, 'r') as f:
-        results_dict = json.load(f)
-
-    message = gen_message(pr,
-                          status['main_stat'],
-                          status['test_stat'],
-                          action_url)
-
-    # Post your comment
-    post_comment(pr, message, match_string='Linter Bot Results:')

From 237f58d1d74481f9e0d6c8e6f7cad3a71a7e5ff2 Mon Sep 17 00:00:00 2001
From: Rocco Meli <r.meli@bluemail.ch>
Date: Wed, 2 Oct 2024 23:20:49 +0200
Subject: [PATCH 2/3] format tables and due

---
 package/MDAnalysis/due.py             | 21 ++++++++-----
 package/MDAnalysis/topology/tables.py | 44 ++++++++++++++++-----------
 package/pyproject.toml                |  7 ++++-
 3 files changed, 46 insertions(+), 26 deletions(-)

diff --git a/package/MDAnalysis/due.py b/package/MDAnalysis/due.py
index 0528bb83e15..7ab1e8b55a9 100644
--- a/package/MDAnalysis/due.py
+++ b/package/MDAnalysis/due.py
@@ -26,30 +26,33 @@
 
 """
 
-__version__ = '0.0.5'
+__version__ = "0.0.5"
 
 
 class InactiveDueCreditCollector(object):
     """Just a stub at the Collector which would not do anything"""
+
     def _donothing(self, *args, **kwargs):
         """Perform no good and no bad"""
-        pass # pylint: disable=unnecessary-pass
+        pass  # pylint: disable=unnecessary-pass
 
     def dcite(self, *args, **kwargs):
         """If I could cite I would"""
+
         def nondecorating_decorator(func):
             return func
+
         return nondecorating_decorator
 
     cite = load = add = _donothing
 
     def __repr__(self):
-        return self.__class__.__name__ + '()'
+        return self.__class__.__name__ + "()"
 
 
 def _donothing_func(*args, **kwargs):
     """Perform no good and no bad"""
-    pass # pylint: disable=unnecessary-pass
+    pass  # pylint: disable=unnecessary-pass
 
 
 try:
@@ -60,17 +63,21 @@ def _donothing_func(*args, **kwargs):
     import duecredit
 
     from duecredit import due, BibTeX, Doi, Url
-    if 'due' in locals() and not hasattr(due, 'cite'):
+
+    if "due" in locals() and not hasattr(due, "cite"):
         raise RuntimeError(
-            "Imported due lacks .cite. DueCredit is now disabled")
+            "Imported due lacks .cite. DueCredit is now disabled"
+        )
 except Exception as err:
     if not isinstance(err, ImportError):
         import logging
         import warnings
+
         errmsg = "Failed to import duecredit due to {}".format(str(err))
         warnings.warn(errmsg)
         logging.getLogger("duecredit").error(
-            "Failed to import duecredit due to {}".format(str(err)))
+            "Failed to import duecredit due to {}".format(str(err))
+        )
     # else:
     #   Do not issue any warnings if duecredit is not installed;
     #   this is the user's choice (Issue #1872)
diff --git a/package/MDAnalysis/topology/tables.py b/package/MDAnalysis/topology/tables.py
index 1c0b51b5817..ea46421cc42 100644
--- a/package/MDAnalysis/topology/tables.py
+++ b/package/MDAnalysis/topology/tables.py
@@ -46,8 +46,10 @@
 .. autodata:: TABLE_VDWRADII
 """
 
+from typing import Any
 
-def kv2dict(s, convertor=str):
+
+def kv2dict(s, convertor: Any = str):
     """Primitive ad-hoc parser of a key-value record list.
 
     * The string *s* should contain each key-value pair on a separate
@@ -172,10 +174,12 @@ def kv2dict(s, convertor=str):
 #: with :func:`MDAnalysis.topology.core.guess_atom_type`.
 atomelements = kv2dict(TABLE_ATOMELEMENTS)
 
+# fmt: off
 elements = ['H',
             'LI', 'BE', 'B', 'C', 'N', 'O', 'F',
             'NA', 'MG', 'AL', 'P', 'SI', 'S', 'CL',
             'K']
+# fmt: on
 
 #: Plain-text table with atomic masses in u.
 TABLE_MASSES = """
@@ -376,28 +380,31 @@ def kv2dict(s, convertor=str):
 #: .. SeeAlso:: :func:`MDAnalysis.topology.core.guess_bonds`
 vdwradii = kv2dict(TABLE_VDWRADII, convertor=float)
 
-Z2SYMB = {1: 'H',                                                                2: 'He', 
-          3: 'Li',   4: 'Be',  5: 'B',   6: 'C',   7: 'N',   8: 'O',   9: 'F',  10: 'Ne', 
-          11: 'Na', 12: 'Mg', 13: 'Al', 14: 'Si', 15: 'P',  16: 'S',  17: 'Cl', 18: 'Ar', 
-          19: 'K',  20: 'Ca', 21: 'Sc', 22: 'Ti', 23: 'V',  24: 'Cr', 25: 'Mn', 26: 'Fe', 
-          27: 'Co', 28: 'Ni', 29: 'Cu', 30: 'Zn', 31: 'Ga', 32: 'Ge', 33: 'As', 34: 'Se', 
-          35: 'Br', 36: 'Kr', 37: 'Rb', 38: 'Sr', 39: 'Y',  40: 'Zr', 41: 'Nb', 42: 'Mo', 
-          43: 'Tc', 44: 'Ru', 45: 'Rh', 46: 'Pd', 47: 'Ag', 48: 'Cd', 49: 'In', 50: 'Sn', 
-          51: 'Sb', 52: 'Te', 53: 'I',  54: 'Xe', 55: 'Cs', 56: 'Ba', 57: 'La', 58: 'Ce', 
-          59: 'Pr', 60: 'Nd', 61: 'Pm', 62: 'Sm', 63: 'Eu', 64: 'Gd', 65: 'Tb', 66: 'Dy', 
-          67: 'Ho', 68: 'Er', 69: 'Tm', 70: 'Yb', 71: 'Lu', 72: 'Hf', 73: 'Ta', 74: 'W', 
-          75: 'Re', 76: 'Os', 77: 'Ir', 78: 'Pt', 79: 'Au', 80: 'Hg', 81: 'Tl', 82: 'Pb', 
-          83: 'Bi', 84: 'Po', 85: 'At', 86: 'Rn', 87: 'Fr', 88: 'Ra', 89: 'Ac', 90: 'Th', 
-          91: 'Pa', 92: 'U',  93: 'Np', 94: 'Pu', 95: 'Am', 96: 'Cm', 97: 'Bk', 98: 'Cf', 
-          99: 'Es', 100: 'Fm', 101: 'Md', 102: 'No', 103: 'Lr', 104: 'Rf', 105: 'Db', 
-          106: 'Sg', 107: 'Bh', 108: 'Hs', 109: 'Mt', 110: 'Ds', 111: 'Rg', 112: 'Cn', 
+# fmt: off
+Z2SYMB = {1: 'H',                                                                2: 'He',
+          3: 'Li',   4: 'Be',  5: 'B',   6: 'C',   7: 'N',   8: 'O',   9: 'F',  10: 'Ne',
+          11: 'Na', 12: 'Mg', 13: 'Al', 14: 'Si', 15: 'P',  16: 'S',  17: 'Cl', 18: 'Ar',
+          19: 'K',  20: 'Ca', 21: 'Sc', 22: 'Ti', 23: 'V',  24: 'Cr', 25: 'Mn', 26: 'Fe',
+          27: 'Co', 28: 'Ni', 29: 'Cu', 30: 'Zn', 31: 'Ga', 32: 'Ge', 33: 'As', 34: 'Se',
+          35: 'Br', 36: 'Kr', 37: 'Rb', 38: 'Sr', 39: 'Y',  40: 'Zr', 41: 'Nb', 42: 'Mo',
+          43: 'Tc', 44: 'Ru', 45: 'Rh', 46: 'Pd', 47: 'Ag', 48: 'Cd', 49: 'In', 50: 'Sn',
+          51: 'Sb', 52: 'Te', 53: 'I',  54: 'Xe', 55: 'Cs', 56: 'Ba', 57: 'La', 58: 'Ce',
+          59: 'Pr', 60: 'Nd', 61: 'Pm', 62: 'Sm', 63: 'Eu', 64: 'Gd', 65: 'Tb', 66: 'Dy',
+          67: 'Ho', 68: 'Er', 69: 'Tm', 70: 'Yb', 71: 'Lu', 72: 'Hf', 73: 'Ta', 74: 'W',
+          75: 'Re', 76: 'Os', 77: 'Ir', 78: 'Pt', 79: 'Au', 80: 'Hg', 81: 'Tl', 82: 'Pb',
+          83: 'Bi', 84: 'Po', 85: 'At', 86: 'Rn', 87: 'Fr', 88: 'Ra', 89: 'Ac', 90: 'Th',
+          91: 'Pa', 92: 'U',  93: 'Np', 94: 'Pu', 95: 'Am', 96: 'Cm', 97: 'Bk', 98: 'Cf',
+          99: 'Es', 100: 'Fm', 101: 'Md', 102: 'No', 103: 'Lr', 104: 'Rf', 105: 'Db',
+          106: 'Sg', 107: 'Bh', 108: 'Hs', 109: 'Mt', 110: 'Ds', 111: 'Rg', 112: 'Cn',
           113: 'Nh', 114: 'Fl', 115: 'Mc', 116: 'Lv', 117: 'Ts', 118: 'Og'}
+# fmt: on
 
-SYMB2Z = {v:k for k, v in Z2SYMB.items()}
+SYMB2Z = {v: k for k, v in Z2SYMB.items()}
 
 # Conversion between SYBYL atom types and corresponding elements
 # Tripos MOL2 file format:
 #   https://web.archive.org/web/*/http://chemyang.ccnu.edu.cn/ccb/server/AIMMS/mol2.pdf
+# fmt: off
 SYBYL2SYMB = {
     "H": "H", "H.spc": "H", "H.t3p": "H",
     "C.3": "C", "C.2": "C", "C.1": "C", "C.ar": "C", "C.cat": "C",
@@ -428,4 +435,5 @@ def kv2dict(s, convertor=str):
     "Se": "Se",
     "Mo": "Mo",
     "Sn": "Sn",
-}
\ No newline at end of file
+}
+# fmt: on
diff --git a/package/pyproject.toml b/package/pyproject.toml
index 3d1ea04bddb..297ed6f3040 100644
--- a/package/pyproject.toml
+++ b/package/pyproject.toml
@@ -125,5 +125,10 @@ MDAnalysis = [
 [tool.black]
 line-length = 79
 target-version = ['py310', 'py311', 'py312']
-extend-exclude = '.'
+include = '''
+(
+tables\.py
+| due\.py
+)
+'''
 required-version = '24'

From 19c3afcf11cec32f5470bbfd95cea00ae9ec377e Mon Sep 17 00:00:00 2001
From: Rocco Meli <r.meli@bluemail.ch>
Date: Thu, 3 Oct 2024 15:20:52 +0200
Subject: [PATCH 3/3] setup.py and visualization

---
 package/MDAnalysis/visualization/__init__.py  |   2 +-
 .../MDAnalysis/visualization/streamlines.py   | 253 ++++++---
 .../visualization/streamlines_3D.py           | 446 ++++++++++-----
 package/pyproject.toml                        |   3 +
 package/setup.py                              | 522 ++++++++++--------
 5 files changed, 799 insertions(+), 427 deletions(-)

diff --git a/package/MDAnalysis/visualization/__init__.py b/package/MDAnalysis/visualization/__init__.py
index 72e488329a2..e113eb73d7b 100644
--- a/package/MDAnalysis/visualization/__init__.py
+++ b/package/MDAnalysis/visualization/__init__.py
@@ -24,4 +24,4 @@
 from . import streamlines
 from . import streamlines_3D
 
-__all__ = ['streamlines', 'streamlines_3D']
+__all__ = ["streamlines", "streamlines_3D"]
diff --git a/package/MDAnalysis/visualization/streamlines.py b/package/MDAnalysis/visualization/streamlines.py
index 965074a43d7..ea9a2bb6342 100644
--- a/package/MDAnalysis/visualization/streamlines.py
+++ b/package/MDAnalysis/visualization/streamlines.py
@@ -52,16 +52,15 @@
     import matplotlib.path
 except ImportError:
     raise ImportError(
-            '2d streamplot module requires: matplotlib.path for its '
-            'path.Path.contains_points method. The installation '
-            'instructions for the matplotlib module can be found here: '
-            'http://matplotlib.org/faq/installing_faq.html?highlight=install'
-            ) from None
+        "2d streamplot module requires: matplotlib.path for its "
+        "path.Path.contains_points method. The installation "
+        "instructions for the matplotlib module can be found here: "
+        "http://matplotlib.org/faq/installing_faq.html?highlight=install"
+    ) from None
 
 import MDAnalysis
 
 
-
 def produce_grid(tuple_of_limits, grid_spacing):
     """Produce a 2D grid for the simulation system.
 
@@ -120,12 +119,16 @@ def split_grid(grid, num_cores):
     # produce an array containing the cartesian coordinates of all vertices in the grid:
     x_array, y_array = grid
     grid_vertex_cartesian_array = np.dstack((x_array, y_array))
-    #the grid_vertex_cartesian_array has N_rows, with each row corresponding to a column of coordinates in the grid (
+    # the grid_vertex_cartesian_array has N_rows, with each row corresponding to a column of coordinates in the grid (
     # so a given row has shape N_rows, 2); overall shape (N_columns_in_grid, N_rows_in_a_column, 2)
-    #although I'll eventually want a pure numpy/scipy/vector-based solution, for now I'll allow loops to simplify the
+    # although I'll eventually want a pure numpy/scipy/vector-based solution, for now I'll allow loops to simplify the
     #  division of the cartesian coordinates into a list of the squares in the grid
-    list_all_squares_in_grid = []  # should eventually be a nested list of all the square vertices in the grid/system
-    list_parent_index_values = []  # want an ordered list of assignment indices for reconstructing the grid positions
+    list_all_squares_in_grid = (
+        []
+    )  # should eventually be a nested list of all the square vertices in the grid/system
+    list_parent_index_values = (
+        []
+    )  # want an ordered list of assignment indices for reconstructing the grid positions
     # in the parent process
     current_column = 0
     while current_column < grid_vertex_cartesian_array.shape[0] - 1:
@@ -134,100 +137,182 @@ def split_grid(grid, num_cores):
         current_row = 0
         while current_row < grid_vertex_cartesian_array.shape[1] - 1:
             # all rows except the top row, which doesn't have a row above it for forming squares
-            bottom_left_vertex_current_square = grid_vertex_cartesian_array[current_column, current_row]
-            bottom_right_vertex_current_square = grid_vertex_cartesian_array[current_column + 1, current_row]
-            top_right_vertex_current_square = grid_vertex_cartesian_array[current_column + 1, current_row + 1]
-            top_left_vertex_current_square = grid_vertex_cartesian_array[current_column, current_row + 1]
-            #append the vertices of this square to the overall list of square vertices:
+            bottom_left_vertex_current_square = grid_vertex_cartesian_array[
+                current_column, current_row
+            ]
+            bottom_right_vertex_current_square = grid_vertex_cartesian_array[
+                current_column + 1, current_row
+            ]
+            top_right_vertex_current_square = grid_vertex_cartesian_array[
+                current_column + 1, current_row + 1
+            ]
+            top_left_vertex_current_square = grid_vertex_cartesian_array[
+                current_column, current_row + 1
+            ]
+            # append the vertices of this square to the overall list of square vertices:
             list_all_squares_in_grid.append(
-                [bottom_left_vertex_current_square, bottom_right_vertex_current_square, top_right_vertex_current_square,
-                    top_left_vertex_current_square])
+                [
+                    bottom_left_vertex_current_square,
+                    bottom_right_vertex_current_square,
+                    top_right_vertex_current_square,
+                    top_left_vertex_current_square,
+                ]
+            )
             list_parent_index_values.append([current_row, current_column])
             current_row += 1
         current_column += 1
-    #split the list of square vertices [[v1,v2,v3,v4],[v1,v2,v3,v4],...,...] into roughly equally-sized sublists to
+    # split the list of square vertices [[v1,v2,v3,v4],[v1,v2,v3,v4],...,...] into roughly equally-sized sublists to
     # be distributed over the available cores on the system:
-    list_square_vertex_arrays_per_core = np.array_split(list_all_squares_in_grid, num_cores)
-    list_parent_index_values = np.array_split(list_parent_index_values, num_cores)
-    return [list_square_vertex_arrays_per_core, list_parent_index_values, current_row, current_column]
-
-
-def per_core_work(topology_file_path, trajectory_file_path, list_square_vertex_arrays_this_core, MDA_selection,
-                  start_frame, end_frame, reconstruction_index_list, maximum_delta_magnitude):
+    list_square_vertex_arrays_per_core = np.array_split(
+        list_all_squares_in_grid, num_cores
+    )
+    list_parent_index_values = np.array_split(
+        list_parent_index_values, num_cores
+    )
+    return [
+        list_square_vertex_arrays_per_core,
+        list_parent_index_values,
+        current_row,
+        current_column,
+    ]
+
+
+def per_core_work(
+    topology_file_path,
+    trajectory_file_path,
+    list_square_vertex_arrays_this_core,
+    MDA_selection,
+    start_frame,
+    end_frame,
+    reconstruction_index_list,
+    maximum_delta_magnitude,
+):
     """Run the analysis on one core.
 
     The code to perform on a given core given the list of square vertices assigned to it.
     """
     # obtain the relevant coordinates for particles of interest
-    universe_object = MDAnalysis.Universe(topology_file_path, trajectory_file_path)
+    universe_object = MDAnalysis.Universe(
+        topology_file_path, trajectory_file_path
+    )
     list_previous_frame_centroids = []
     list_previous_frame_indices = []
-    #define some utility functions for trajectory iteration:
+    # define some utility functions for trajectory iteration:
 
     def produce_list_indices_point_in_polygon_this_frame(vertex_coord_list):
         list_indices_point_in_polygon = []
         for square_vertices in vertex_coord_list:
             path_object = matplotlib.path.Path(square_vertices)
-            index_list_in_polygon = np.where(path_object.contains_points(relevant_particle_coordinate_array_xy))
+            index_list_in_polygon = np.where(
+                path_object.contains_points(
+                    relevant_particle_coordinate_array_xy
+                )
+            )
             list_indices_point_in_polygon.append(index_list_in_polygon)
         return list_indices_point_in_polygon
 
     def produce_list_centroids_this_frame(list_indices_in_polygon):
         list_centroids_this_frame = []
         for indices in list_indices_in_polygon:
-            if not indices[0].size > 0:  # if there are no particles of interest in this particular square
+            if (
+                not indices[0].size > 0
+            ):  # if there are no particles of interest in this particular square
                 list_centroids_this_frame.append(None)
             else:
-                current_coordinate_array_in_square = relevant_particle_coordinate_array_xy[indices]
-                current_square_indices_centroid = np.average(current_coordinate_array_in_square, axis=0)
-                list_centroids_this_frame.append(current_square_indices_centroid)
+                current_coordinate_array_in_square = (
+                    relevant_particle_coordinate_array_xy[indices]
+                )
+                current_square_indices_centroid = np.average(
+                    current_coordinate_array_in_square, axis=0
+                )
+                list_centroids_this_frame.append(
+                    current_square_indices_centroid
+                )
         return list_centroids_this_frame  # a list of numpy xy centroid arrays for this frame
 
     for ts in universe_object.trajectory:
         if ts.frame < start_frame:  # don't start until first specified frame
             continue
-        relevant_particle_coordinate_array_xy = universe_object.select_atoms(MDA_selection).positions[..., :-1]
+        relevant_particle_coordinate_array_xy = universe_object.select_atoms(
+            MDA_selection
+        ).positions[..., :-1]
         # only 2D / xy coords for now
-        #I will need a list of indices for relevant particles falling within each square in THIS frame:
-        list_indices_in_squares_this_frame = produce_list_indices_point_in_polygon_this_frame(
-            list_square_vertex_arrays_this_core)
-        #likewise, I will need a list of centroids of particles in each square (same order as above list):
-        list_centroids_in_squares_this_frame = produce_list_centroids_this_frame(list_indices_in_squares_this_frame)
-        if list_previous_frame_indices:  # if the previous frame had indices in at least one square I will need to use
+        # I will need a list of indices for relevant particles falling within each square in THIS frame:
+        list_indices_in_squares_this_frame = (
+            produce_list_indices_point_in_polygon_this_frame(
+                list_square_vertex_arrays_this_core
+            )
+        )
+        # likewise, I will need a list of centroids of particles in each square (same order as above list):
+        list_centroids_in_squares_this_frame = (
+            produce_list_centroids_this_frame(
+                list_indices_in_squares_this_frame
+            )
+        )
+        if (
+            list_previous_frame_indices
+        ):  # if the previous frame had indices in at least one square I will need to use
             #  those indices to generate the updates to the corresponding centroids in this frame:
-            list_centroids_this_frame_using_indices_from_last_frame = produce_list_centroids_this_frame(
-                list_previous_frame_indices)
-            #I need to write a velocity of zero if there are any 'empty' squares in either frame:
+            list_centroids_this_frame_using_indices_from_last_frame = (
+                produce_list_centroids_this_frame(list_previous_frame_indices)
+            )
+            # I need to write a velocity of zero if there are any 'empty' squares in either frame:
             xy_deltas_to_write = []
-            for square_1_centroid, square_2_centroid in zip(list_centroids_this_frame_using_indices_from_last_frame,
-                                                            list_previous_frame_centroids):
+            for square_1_centroid, square_2_centroid in zip(
+                list_centroids_this_frame_using_indices_from_last_frame,
+                list_previous_frame_centroids,
+            ):
                 if square_1_centroid is None or square_2_centroid is None:
                     xy_deltas_to_write.append([0, 0])
                 else:
-                    xy_deltas_to_write.append(np.subtract(square_1_centroid, square_2_centroid).tolist())
+                    xy_deltas_to_write.append(
+                        np.subtract(
+                            square_1_centroid, square_2_centroid
+                        ).tolist()
+                    )
 
-            #xy_deltas_to_write = np.subtract(np.array(
+            # xy_deltas_to_write = np.subtract(np.array(
             # list_centroids_this_frame_using_indices_from_last_frame),np.array(list_previous_frame_centroids))
             xy_deltas_to_write = np.array(xy_deltas_to_write)
-            #now filter the array to only contain distances in the range [-8,8] as a placeholder for dealing with PBC
+            # now filter the array to only contain distances in the range [-8,8] as a placeholder for dealing with PBC
             #  issues (Matthieu seemed to use a limit of 8 as well);
-            xy_deltas_to_write = np.clip(xy_deltas_to_write, -maximum_delta_magnitude, maximum_delta_magnitude)
+            xy_deltas_to_write = np.clip(
+                xy_deltas_to_write,
+                -maximum_delta_magnitude,
+                maximum_delta_magnitude,
+            )
 
-            #with the xy and dx,dy values calculated I need to set the values from this frame to previous frame
+            # with the xy and dx,dy values calculated I need to set the values from this frame to previous frame
             # values in anticipation of the next frame:
-            list_previous_frame_centroids = list_centroids_in_squares_this_frame[:]
+            list_previous_frame_centroids = (
+                list_centroids_in_squares_this_frame[:]
+            )
             list_previous_frame_indices = list_indices_in_squares_this_frame[:]
         else:  # either no points in squares or after the first frame I'll just reset the 'previous' values so they
             # can be used when consecutive frames have proper values
-            list_previous_frame_centroids = list_centroids_in_squares_this_frame[:]
+            list_previous_frame_centroids = (
+                list_centroids_in_squares_this_frame[:]
+            )
             list_previous_frame_indices = list_indices_in_squares_this_frame[:]
         if ts.frame > end_frame:
             break  # stop here
     return list(zip(reconstruction_index_list, xy_deltas_to_write.tolist()))
 
 
-def generate_streamlines(topology_file_path, trajectory_file_path, grid_spacing, MDA_selection, start_frame,
-                         end_frame, xmin, xmax, ymin, ymax, maximum_delta_magnitude, num_cores='maximum'):
+def generate_streamlines(
+    topology_file_path,
+    trajectory_file_path,
+    grid_spacing,
+    MDA_selection,
+    start_frame,
+    end_frame,
+    xmin,
+    xmax,
+    ymin,
+    ymax,
+    maximum_delta_magnitude,
+    num_cores="maximum",
+):
     r"""Produce the x and y components of a 2D streamplot data set.
 
     Parameters
@@ -311,35 +396,58 @@ def generate_streamlines(topology_file_path, trajectory_file_path, grid_spacing,
 
     """
     # work out the number of cores to use:
-    if num_cores == 'maximum':
+    if num_cores == "maximum":
         num_cores = multiprocessing.cpu_count()  # use all available cores
     else:
         num_cores = num_cores  # use the value specified by the user
-        #assert isinstance(num_cores,(int,long)), "The number of specified cores must (of course) be an integer."
-    np.seterr(all='warn', over='raise')
+        # assert isinstance(num_cores,(int,long)), "The number of specified cores must (of course) be an integer."
+    np.seterr(all="warn", over="raise")
     parent_list_deltas = []  # collect all data from child processes here
 
     def log_result_to_parent(delta_array):
         parent_list_deltas.extend(delta_array)
 
     tuple_of_limits = (xmin, xmax, ymin, ymax)
-    grid = produce_grid(tuple_of_limits=tuple_of_limits, grid_spacing=grid_spacing)
-    list_square_vertex_arrays_per_core, list_parent_index_values, total_rows, total_columns = \
-        split_grid(grid=grid,
-                   num_cores=num_cores)
+    grid = produce_grid(
+        tuple_of_limits=tuple_of_limits, grid_spacing=grid_spacing
+    )
+    (
+        list_square_vertex_arrays_per_core,
+        list_parent_index_values,
+        total_rows,
+        total_columns,
+    ) = split_grid(grid=grid, num_cores=num_cores)
     pool = multiprocessing.Pool(num_cores)
-    for vertex_sublist, index_sublist in zip(list_square_vertex_arrays_per_core, list_parent_index_values):
-        pool.apply_async(per_core_work, args=(
-            topology_file_path, trajectory_file_path, vertex_sublist, MDA_selection, start_frame, end_frame,
-            index_sublist, maximum_delta_magnitude), callback=log_result_to_parent)
+    for vertex_sublist, index_sublist in zip(
+        list_square_vertex_arrays_per_core, list_parent_index_values
+    ):
+        pool.apply_async(
+            per_core_work,
+            args=(
+                topology_file_path,
+                trajectory_file_path,
+                vertex_sublist,
+                MDA_selection,
+                start_frame,
+                end_frame,
+                index_sublist,
+                maximum_delta_magnitude,
+            ),
+            callback=log_result_to_parent,
+        )
     pool.close()
     pool.join()
     dx_array = np.zeros((total_rows, total_columns))
     dy_array = np.zeros((total_rows, total_columns))
-    #the parent_list_deltas is shaped like this: [ ([row_index,column_index],[dx,dy]), ... (...),...,]
-    for index_array, delta_array in parent_list_deltas:  # go through the list in the parent process and assign to the
+    # the parent_list_deltas is shaped like this: [ ([row_index,column_index],[dx,dy]), ... (...),...,]
+    for (
+        index_array,
+        delta_array,
+    ) in (
+        parent_list_deltas
+    ):  # go through the list in the parent process and assign to the
         #  appropriate positions in the dx and dy matrices:
-        #build in a filter to replace all values at the cap (currently between -8,8) with 0 to match Matthieu's code
+        # build in a filter to replace all values at the cap (currently between -8,8) with 0 to match Matthieu's code
         # (I think eventually we'll reduce the cap to a narrower boundary though)
         index_1 = index_array.tolist()[0]
         index_2 = index_array.tolist()[1]
@@ -352,9 +460,14 @@ def log_result_to_parent(delta_array):
         else:
             dy_array[index_1, index_2] = delta_array[1]
 
-    #at Matthieu's request, we now want to calculate the average and standard deviation of the displacement values:
-    displacement_array = np.sqrt(dx_array ** 2 + dy_array ** 2)
+    # at Matthieu's request, we now want to calculate the average and standard deviation of the displacement values:
+    displacement_array = np.sqrt(dx_array**2 + dy_array**2)
     average_displacement = np.average(displacement_array)
     standard_deviation_of_displacement = np.std(displacement_array)
 
-    return (dx_array, dy_array, average_displacement, standard_deviation_of_displacement)
+    return (
+        dx_array,
+        dy_array,
+        average_displacement,
+        standard_deviation_of_displacement,
+    )
diff --git a/package/MDAnalysis/visualization/streamlines_3D.py b/package/MDAnalysis/visualization/streamlines_3D.py
index 1f85851c16a..5c9a03c4e12 100644
--- a/package/MDAnalysis/visualization/streamlines_3D.py
+++ b/package/MDAnalysis/visualization/streamlines_3D.py
@@ -56,8 +56,9 @@
 import MDAnalysis
 
 
-def determine_container_limits(topology_file_path, trajectory_file_path,
-                               buffer_value):
+def determine_container_limits(
+    topology_file_path, trajectory_file_path, buffer_value
+):
     """Calculate the extent of the atom coordinates + buffer.
 
     A function for the parent process which should take the input trajectory
@@ -73,19 +74,29 @@ def determine_container_limits(topology_file_path, trajectory_file_path,
     buffer_value : float
         buffer value (padding) in +/- {x, y, z}
     """
-    universe_object = MDAnalysis.Universe(topology_file_path, trajectory_file_path)
-    all_atom_selection = universe_object.select_atoms('all')  # select all particles
+    universe_object = MDAnalysis.Universe(
+        topology_file_path, trajectory_file_path
+    )
+    all_atom_selection = universe_object.select_atoms(
+        "all"
+    )  # select all particles
     all_atom_coordinate_array = all_atom_selection.positions
     x_min, x_max, y_min, y_max, z_min, z_max = [
         all_atom_coordinate_array[..., 0].min(),
-        all_atom_coordinate_array[..., 0].max(), all_atom_coordinate_array[..., 1].min(),
-        all_atom_coordinate_array[..., 1].max(), all_atom_coordinate_array[..., 2].min(),
-        all_atom_coordinate_array[..., 2].max()]
-    tuple_of_limits = \
-        (
-            x_min - buffer_value,
-            x_max + buffer_value, y_min - buffer_value, y_max + buffer_value, z_min - buffer_value,
-            z_max + buffer_value)  # using buffer_value to catch particles near edges
+        all_atom_coordinate_array[..., 0].max(),
+        all_atom_coordinate_array[..., 1].min(),
+        all_atom_coordinate_array[..., 1].max(),
+        all_atom_coordinate_array[..., 2].min(),
+        all_atom_coordinate_array[..., 2].max(),
+    ]
+    tuple_of_limits = (
+        x_min - buffer_value,
+        x_max + buffer_value,
+        y_min - buffer_value,
+        y_max + buffer_value,
+        z_min - buffer_value,
+        z_max + buffer_value,
+    )  # using buffer_value to catch particles near edges
     return tuple_of_limits
 
 
@@ -109,7 +120,11 @@ def produce_grid(tuple_of_limits, grid_spacing):
 
     """
     x_min, x_max, y_min, y_max, z_min, z_max = tuple_of_limits
-    grid = np.mgrid[x_min:x_max:grid_spacing, y_min:y_max:grid_spacing, z_min:z_max:grid_spacing]
+    grid = np.mgrid[
+        x_min:x_max:grid_spacing,
+        y_min:y_max:grid_spacing,
+        z_min:z_max:grid_spacing,
+    ]
     return grid
 
 
@@ -139,78 +154,124 @@ def split_grid(grid, num_cores):
     num_z_values = z.shape[-1]
     num_sheets = z.shape[0]
     delta_array_shape = tuple(
-        [n - 1 for n in x.shape])  # the final target shape for return delta arrays is n-1 in each dimension
+        [n - 1 for n in x.shape]
+    )  # the final target shape for return delta arrays is n-1 in each dimension
 
     ordered_list_per_sheet_x_values = []
-    for x_sheet in x:  # each x_sheet should have shape (25,23) and the same x value in each element
+    for (
+        x_sheet
+    ) in (
+        x
+    ):  # each x_sheet should have shape (25,23) and the same x value in each element
         array_all_x_values_current_sheet = x_sheet.flatten()
-        ordered_list_per_sheet_x_values.append(array_all_x_values_current_sheet)
+        ordered_list_per_sheet_x_values.append(
+            array_all_x_values_current_sheet
+        )
     ordered_list_per_sheet_y_values = []
     for y_columns in y:
         array_all_y_values_current_sheet = y_columns.flatten()
-        ordered_list_per_sheet_y_values.append(array_all_y_values_current_sheet)
+        ordered_list_per_sheet_y_values.append(
+            array_all_y_values_current_sheet
+        )
     ordered_list_per_sheet_z_values = []
     for z_slices in z:
         array_all_z_values_current_sheet = z_slices.flatten()
-        ordered_list_per_sheet_z_values.append(array_all_z_values_current_sheet)
+        ordered_list_per_sheet_z_values.append(
+            array_all_z_values_current_sheet
+        )
 
     ordered_list_cartesian_coordinates_per_sheet = []
-    for x_sheet_coords, y_sheet_coords, z_sheet_coords in zip(ordered_list_per_sheet_x_values,
-                                                              ordered_list_per_sheet_y_values,
-                                                              ordered_list_per_sheet_z_values):
-        ordered_list_cartesian_coordinates_per_sheet.append(list(zip(x_sheet_coords, y_sheet_coords, z_sheet_coords)))
-    array_ordered_cartesian_coords_per_sheet = np.array(ordered_list_cartesian_coordinates_per_sheet)
-    #now I'm going to want to build cubes in an ordered fashion, and in such a way that I can track the index /
+    for x_sheet_coords, y_sheet_coords, z_sheet_coords in zip(
+        ordered_list_per_sheet_x_values,
+        ordered_list_per_sheet_y_values,
+        ordered_list_per_sheet_z_values,
+    ):
+        ordered_list_cartesian_coordinates_per_sheet.append(
+            list(zip(x_sheet_coords, y_sheet_coords, z_sheet_coords))
+        )
+    array_ordered_cartesian_coords_per_sheet = np.array(
+        ordered_list_cartesian_coordinates_per_sheet
+    )
+    # now I'm going to want to build cubes in an ordered fashion, and in such a way that I can track the index /
     # centroid of each cube for domain decomposition / reconstruction and mayavi mlab.flow() input
-    #cubes will be formed from N - 1 base sheets combined with subsequent sheets
+    # cubes will be formed from N - 1 base sheets combined with subsequent sheets
     current_base_sheet = 0
     dictionary_cubes_centroids_indices = {}
     cube_counter = 0
     while current_base_sheet < num_sheets - 1:
-        current_base_sheet_array = array_ordered_cartesian_coords_per_sheet[current_base_sheet]
+        current_base_sheet_array = array_ordered_cartesian_coords_per_sheet[
+            current_base_sheet
+        ]
         current_top_sheet_array = array_ordered_cartesian_coords_per_sheet[
-            current_base_sheet + 1]  # the points of the sheet 'to the right' in the grid
+            current_base_sheet + 1
+        ]  # the points of the sheet 'to the right' in the grid
         current_index = 0
         while current_index < current_base_sheet_array.shape[0] - num_z_values:
             # iterate through all the indices in each of the sheet arrays (careful to avoid extra
             # points not needed for cubes)
-            column_z_level = 0  # start at the bottom of a given 4-point column and work up
+            column_z_level = (
+                0  # start at the bottom of a given 4-point column and work up
+            )
             while column_z_level < num_z_values - 1:
                 current_list_cube_vertices = []
-                first_two_vertices_base_sheet = current_base_sheet_array[current_index:current_index + 2, ...].tolist()
-                first_two_vertices_top_sheet = current_top_sheet_array[current_index:current_index + 2, ...].tolist()
-                next_two_vertices_base_sheet = current_base_sheet_array[current_index +
-                                                                        num_z_values: 2 +
-                                                                        num_z_values + current_index, ...].tolist()
-                next_two_vertices_top_sheet = current_top_sheet_array[current_index +
-                                                                      num_z_values: 2 +
-                                                                      num_z_values + current_index, ...].tolist()
+                first_two_vertices_base_sheet = current_base_sheet_array[
+                    current_index : current_index + 2, ...
+                ].tolist()
+                first_two_vertices_top_sheet = current_top_sheet_array[
+                    current_index : current_index + 2, ...
+                ].tolist()
+                next_two_vertices_base_sheet = current_base_sheet_array[
+                    current_index
+                    + num_z_values : 2
+                    + num_z_values
+                    + current_index,
+                    ...,
+                ].tolist()
+                next_two_vertices_top_sheet = current_top_sheet_array[
+                    current_index
+                    + num_z_values : 2
+                    + num_z_values
+                    + current_index,
+                    ...,
+                ].tolist()
                 for vertex_set in [
-                    first_two_vertices_base_sheet, first_two_vertices_top_sheet,
-                    next_two_vertices_base_sheet, next_two_vertices_top_sheet
+                    first_two_vertices_base_sheet,
+                    first_two_vertices_top_sheet,
+                    next_two_vertices_base_sheet,
+                    next_two_vertices_top_sheet,
                 ]:
                     current_list_cube_vertices.extend(vertex_set)
                 vertex_array = np.array(current_list_cube_vertices)
-                assert vertex_array.shape == (8, 3), "vertex_array has incorrect shape"
-                cube_centroid = np.average(np.array(current_list_cube_vertices), axis=0)
+                assert vertex_array.shape == (
+                    8,
+                    3,
+                ), "vertex_array has incorrect shape"
+                cube_centroid = np.average(
+                    np.array(current_list_cube_vertices), axis=0
+                )
                 dictionary_cubes_centroids_indices[cube_counter] = {
-                    'centroid': cube_centroid,
-                    'vertex_list': current_list_cube_vertices}
+                    "centroid": cube_centroid,
+                    "vertex_list": current_list_cube_vertices,
+                }
                 cube_counter += 1
                 current_index += 1
                 column_z_level += 1
-                if column_z_level == num_z_values - 1:  # the loop will break but I should also increment the
+                if (
+                    column_z_level == num_z_values - 1
+                ):  # the loop will break but I should also increment the
                     # current_index
                     current_index += 1
         current_base_sheet += 1
     total_cubes = len(dictionary_cubes_centroids_indices)
 
-    #produce an array of pseudo cube indices (actually the dictionary keys which are cube numbers in string format):
+    # produce an array of pseudo cube indices (actually the dictionary keys which are cube numbers in string format):
     pseudo_cube_indices = np.arange(0, total_cubes)
-    sublist_of_cube_indices_per_core = np.array_split(pseudo_cube_indices, num_cores)
-    #now, the split of pseudoindices seems to work well, and the above sublist_of_cube_indices_per_core is a list of
+    sublist_of_cube_indices_per_core = np.array_split(
+        pseudo_cube_indices, num_cores
+    )
+    # now, the split of pseudoindices seems to work well, and the above sublist_of_cube_indices_per_core is a list of
     # arrays of cube numbers / keys in the original dictionary
-    #now I think I'll try to produce a list of dictionaries that each contain their assigned cubes based on the above
+    # now I think I'll try to produce a list of dictionaries that each contain their assigned cubes based on the above
     #  per core split
     list_dictionaries_for_cores = []
     subdictionary_counter = 0
@@ -224,11 +285,22 @@ def split_grid(grid, num_cores):
             items_popped += 1
         list_dictionaries_for_cores.append(current_core_dictionary)
         subdictionary_counter += 1
-    return list_dictionaries_for_cores, total_cubes, num_sheets, delta_array_shape
-
-
-def per_core_work(start_frame_coord_array, end_frame_coord_array, dictionary_cube_data_this_core, MDA_selection,
-                  start_frame, end_frame):
+    return (
+        list_dictionaries_for_cores,
+        total_cubes,
+        num_sheets,
+        delta_array_shape,
+    )
+
+
+def per_core_work(
+    start_frame_coord_array,
+    end_frame_coord_array,
+    dictionary_cube_data_this_core,
+    MDA_selection,
+    start_frame,
+    end_frame,
+):
     """Run the analysis on one core.
 
     The code to perform on a given core given the dictionary of cube data.
@@ -237,82 +309,137 @@ def per_core_work(start_frame_coord_array, end_frame_coord_array, dictionary_cub
     list_previous_frame_indices = []
     # define some utility functions for trajectory iteration:
 
-    def point_in_cube(array_point_coordinates, list_cube_vertices, cube_centroid):
+    def point_in_cube(
+        array_point_coordinates, list_cube_vertices, cube_centroid
+    ):
         """Determine if an array of coordinates are within a cube."""
-        #the simulation particle point can't be more than half the cube side length away from the cube centroid in
+        # the simulation particle point can't be more than half the cube side length away from the cube centroid in
         # any given dimension:
         array_cube_vertices = np.array(list_cube_vertices)
-        cube_half_side_length = scipy.spatial.distance.pdist(array_cube_vertices, 'euclidean').min() / 2.0
-        array_cube_vertex_distances_from_centroid = scipy.spatial.distance.cdist(array_cube_vertices,
-                                                                                 cube_centroid[np.newaxis, :])
-        np.testing.assert_allclose(array_cube_vertex_distances_from_centroid.min(),
-                                          array_cube_vertex_distances_from_centroid.max(), rtol=0, atol=1.5e-4,
-                                          err_msg="not all cube vertex to centroid distances are the same, "
-                                                  "so not a true cube")
-        absolute_delta_coords = np.absolute(np.subtract(array_point_coordinates, cube_centroid))
+        cube_half_side_length = (
+            scipy.spatial.distance.pdist(
+                array_cube_vertices, "euclidean"
+            ).min()
+            / 2.0
+        )
+        array_cube_vertex_distances_from_centroid = (
+            scipy.spatial.distance.cdist(
+                array_cube_vertices, cube_centroid[np.newaxis, :]
+            )
+        )
+        np.testing.assert_allclose(
+            array_cube_vertex_distances_from_centroid.min(),
+            array_cube_vertex_distances_from_centroid.max(),
+            rtol=0,
+            atol=1.5e-4,
+            err_msg="not all cube vertex to centroid distances are the same, "
+            "so not a true cube",
+        )
+        absolute_delta_coords = np.absolute(
+            np.subtract(array_point_coordinates, cube_centroid)
+        )
         absolute_delta_x_coords = absolute_delta_coords[..., 0]
-        indices_delta_x_acceptable = np.where(absolute_delta_x_coords <= cube_half_side_length)
+        indices_delta_x_acceptable = np.where(
+            absolute_delta_x_coords <= cube_half_side_length
+        )
         absolute_delta_y_coords = absolute_delta_coords[..., 1]
-        indices_delta_y_acceptable = np.where(absolute_delta_y_coords <= cube_half_side_length)
+        indices_delta_y_acceptable = np.where(
+            absolute_delta_y_coords <= cube_half_side_length
+        )
         absolute_delta_z_coords = absolute_delta_coords[..., 2]
-        indices_delta_z_acceptable = np.where(absolute_delta_z_coords <= cube_half_side_length)
-        intersection_xy_acceptable_arrays = np.intersect1d(indices_delta_x_acceptable[0],
-                                                              indices_delta_y_acceptable[0])
-        overall_indices_points_in_current_cube = np.intersect1d(intersection_xy_acceptable_arrays,
-                                                                   indices_delta_z_acceptable[0])
+        indices_delta_z_acceptable = np.where(
+            absolute_delta_z_coords <= cube_half_side_length
+        )
+        intersection_xy_acceptable_arrays = np.intersect1d(
+            indices_delta_x_acceptable[0], indices_delta_y_acceptable[0]
+        )
+        overall_indices_points_in_current_cube = np.intersect1d(
+            intersection_xy_acceptable_arrays, indices_delta_z_acceptable[0]
+        )
         return overall_indices_points_in_current_cube
 
-    def update_dictionary_point_in_cube_start_frame(array_simulation_particle_coordinates,
-                                                    dictionary_cube_data_this_core):
+    def update_dictionary_point_in_cube_start_frame(
+        array_simulation_particle_coordinates, dictionary_cube_data_this_core
+    ):
         """Basically update the cube dictionary objects assigned to this core to contain a new key/value pair
         corresponding to the indices of the relevant particles that fall within a given cube. Also, for a given cube,
-        store a key/value pair for the centroid of the particles that fall within the cube."""
+        store a key/value pair for the centroid of the particles that fall within the cube.
+        """
         cube_counter = 0
         for key, cube in dictionary_cube_data_this_core.items():
-            index_list_in_cube = point_in_cube(array_simulation_particle_coordinates, cube['vertex_list'],
-                                               cube['centroid'])
-            cube['start_frame_index_list_in_cube'] = index_list_in_cube
-            if len(index_list_in_cube) > 0:  # if there's at least one particle in this cube
-                centroid_particles_in_cube = np.average(array_simulation_particle_coordinates[index_list_in_cube],
-                                                           axis=0)
-                cube['centroid_of_particles_first_frame'] = centroid_particles_in_cube
+            index_list_in_cube = point_in_cube(
+                array_simulation_particle_coordinates,
+                cube["vertex_list"],
+                cube["centroid"],
+            )
+            cube["start_frame_index_list_in_cube"] = index_list_in_cube
+            if (
+                len(index_list_in_cube) > 0
+            ):  # if there's at least one particle in this cube
+                centroid_particles_in_cube = np.average(
+                    array_simulation_particle_coordinates[index_list_in_cube],
+                    axis=0,
+                )
+                cube["centroid_of_particles_first_frame"] = (
+                    centroid_particles_in_cube
+                )
             else:  # empty cube
-                cube['centroid_of_particles_first_frame'] = None
+                cube["centroid_of_particles_first_frame"] = None
             cube_counter += 1
 
-    def update_dictionary_end_frame(array_simulation_particle_coordinates, dictionary_cube_data_this_core):
+    def update_dictionary_end_frame(
+        array_simulation_particle_coordinates, dictionary_cube_data_this_core
+    ):
         """Update the cube dictionary objects again as appropriate for the second and final frame."""
         cube_counter = 0
         for key, cube in dictionary_cube_data_this_core.items():
             # if there were no particles in the cube in the first frame, then set dx,dy,dz each to 0
-            if cube['centroid_of_particles_first_frame'] is None:
-                cube['dx'] = 0
-                cube['dy'] = 0
-                cube['dz'] = 0
+            if cube["centroid_of_particles_first_frame"] is None:
+                cube["dx"] = 0
+                cube["dy"] = 0
+                cube["dz"] = 0
             else:  # there was at least one particle in the starting cube so we can get dx,dy,dz centroid values
-                new_coordinate_array_for_particles_starting_in_this_cube = array_simulation_particle_coordinates[
-                    cube['start_frame_index_list_in_cube']]
+                new_coordinate_array_for_particles_starting_in_this_cube = (
+                    array_simulation_particle_coordinates[
+                        cube["start_frame_index_list_in_cube"]
+                    ]
+                )
                 new_centroid_for_particles_starting_in_this_cube = np.average(
-                    new_coordinate_array_for_particles_starting_in_this_cube, axis=0)
-                cube['centroid_of_paticles_final_frame'] = new_centroid_for_particles_starting_in_this_cube
-                delta_centroid_array_this_cube = new_centroid_for_particles_starting_in_this_cube - cube[
-                    'centroid_of_particles_first_frame']
-                cube['dx'] = delta_centroid_array_this_cube[0]
-                cube['dy'] = delta_centroid_array_this_cube[1]
-                cube['dz'] = delta_centroid_array_this_cube[2]
+                    new_coordinate_array_for_particles_starting_in_this_cube,
+                    axis=0,
+                )
+                cube["centroid_of_paticles_final_frame"] = (
+                    new_centroid_for_particles_starting_in_this_cube
+                )
+                delta_centroid_array_this_cube = (
+                    new_centroid_for_particles_starting_in_this_cube
+                    - cube["centroid_of_particles_first_frame"]
+                )
+                cube["dx"] = delta_centroid_array_this_cube[0]
+                cube["dy"] = delta_centroid_array_this_cube[1]
+                cube["dz"] = delta_centroid_array_this_cube[2]
             cube_counter += 1
 
-    #now that the parent process is dealing with the universe object & grabbing required coordinates, each child
+    # now that the parent process is dealing with the universe object & grabbing required coordinates, each child
     # process only needs to take the coordinate arrays & perform the operations with its assigned cubes (no more file
     #  opening and trajectory iteration on each core--which I'm hoping will substantially reduce the physical memory
     # footprint of my 3D streamplot code)
-    update_dictionary_point_in_cube_start_frame(start_frame_coord_array, dictionary_cube_data_this_core)
-    update_dictionary_end_frame(end_frame_coord_array, dictionary_cube_data_this_core)
+    update_dictionary_point_in_cube_start_frame(
+        start_frame_coord_array, dictionary_cube_data_this_core
+    )
+    update_dictionary_end_frame(
+        end_frame_coord_array, dictionary_cube_data_this_core
+    )
     return dictionary_cube_data_this_core
 
 
-def produce_coordinate_arrays_single_process(topology_file_path, trajectory_file_path, MDA_selection, start_frame,
-                                             end_frame):
+def produce_coordinate_arrays_single_process(
+    topology_file_path,
+    trajectory_file_path,
+    MDA_selection,
+    start_frame,
+    end_frame,
+):
     """Generate coordinate arrays.
 
     To reduce memory footprint produce only a single MDA selection and get
@@ -321,24 +448,46 @@ def produce_coordinate_arrays_single_process(topology_file_path, trajectory_file
     waste memory.
 
     """
-    universe_object = MDAnalysis.Universe(topology_file_path, trajectory_file_path)
+    universe_object = MDAnalysis.Universe(
+        topology_file_path, trajectory_file_path
+    )
     relevant_particles = universe_object.select_atoms(MDA_selection)
     # pull out coordinate arrays from desired frames:
     for ts in universe_object.trajectory:
         if ts.frame > end_frame:
             break  # stop here
         if ts.frame == start_frame:
-            start_frame_relevant_particle_coordinate_array_xyz = relevant_particles.positions
+            start_frame_relevant_particle_coordinate_array_xyz = (
+                relevant_particles.positions
+            )
         elif ts.frame == end_frame:
-            end_frame_relevant_particle_coordinate_array_xyz = relevant_particles.positions
+            end_frame_relevant_particle_coordinate_array_xyz = (
+                relevant_particles.positions
+            )
         else:
             continue
-    return (start_frame_relevant_particle_coordinate_array_xyz, end_frame_relevant_particle_coordinate_array_xyz)
-
-
-def generate_streamlines_3d(topology_file_path, trajectory_file_path, grid_spacing, MDA_selection, start_frame,
-                            end_frame, xmin, xmax, ymin, ymax, zmin, zmax, maximum_delta_magnitude=2.0,
-                            num_cores='maximum'):
+    return (
+        start_frame_relevant_particle_coordinate_array_xyz,
+        end_frame_relevant_particle_coordinate_array_xyz,
+    )
+
+
+def generate_streamlines_3d(
+    topology_file_path,
+    trajectory_file_path,
+    grid_spacing,
+    MDA_selection,
+    start_frame,
+    end_frame,
+    xmin,
+    xmax,
+    ymin,
+    ymax,
+    zmin,
+    zmax,
+    maximum_delta_magnitude=2.0,
+    num_cores="maximum",
+):
     r"""Produce the x, y and z components of a 3D streamplot data set.
 
     Parameters
@@ -439,68 +588,91 @@ def generate_streamlines_3d(topology_file_path, trajectory_file_path, grid_spaci
     .. _mayavi: http://docs.enthought.com/mayavi/mayavi/
     """
     # work out the number of cores to use:
-    if num_cores == 'maximum':
+    if num_cores == "maximum":
         num_cores = multiprocessing.cpu_count()  # use all available cores
     else:
         num_cores = num_cores  # use the value specified by the user
         # assert isinstance(num_cores,(int,long)), "The number of specified cores must (of course) be an integer."
-    np.seterr(all='warn', over='raise')
+    np.seterr(all="warn", over="raise")
     parent_cube_dictionary = {}  # collect all data from child processes here
 
     def log_result_to_parent(process_dict):
         parent_cube_dictionary.update(process_dict)
 
-    #step 1: produce tuple of cartesian coordinate limits for the first frame
-    #tuple_of_limits = determine_container_limits(topology_file_path = topology_file_path,trajectory_file_path =
+    # step 1: produce tuple of cartesian coordinate limits for the first frame
+    # tuple_of_limits = determine_container_limits(topology_file_path = topology_file_path,trajectory_file_path =
     # trajectory_file_path,buffer_value=buffer_value)
     tuple_of_limits = (xmin, xmax, ymin, ymax, zmin, zmax)
-    #step 2: produce a suitable grid (will assume that grid size / container size does not vary during simulation--or
+    # step 2: produce a suitable grid (will assume that grid size / container size does not vary during simulation--or
     #  at least not beyond the buffer limit, such that this grid can be used for all subsequent frames)
-    grid = produce_grid(tuple_of_limits=tuple_of_limits, grid_spacing=grid_spacing)
-    #step 3: split the grid into a dictionary of cube information that can be sent to each core for processing:
-    list_dictionaries_for_cores, total_cubes, num_sheets, delta_array_shape = split_grid(grid=grid, num_cores=num_cores)
-    #step 3b: produce required coordinate arrays on a single core to avoid making a universe object on each core:
-    start_frame_coord_array, end_frame_coord_array = produce_coordinate_arrays_single_process(topology_file_path,
-                                                                                              trajectory_file_path,
-                                                                                              MDA_selection,
-                                                                                              start_frame, end_frame)
-    #step 4: per process work using the above grid data split
+    grid = produce_grid(
+        tuple_of_limits=tuple_of_limits, grid_spacing=grid_spacing
+    )
+    # step 3: split the grid into a dictionary of cube information that can be sent to each core for processing:
+    list_dictionaries_for_cores, total_cubes, num_sheets, delta_array_shape = (
+        split_grid(grid=grid, num_cores=num_cores)
+    )
+    # step 3b: produce required coordinate arrays on a single core to avoid making a universe object on each core:
+    start_frame_coord_array, end_frame_coord_array = (
+        produce_coordinate_arrays_single_process(
+            topology_file_path,
+            trajectory_file_path,
+            MDA_selection,
+            start_frame,
+            end_frame,
+        )
+    )
+    # step 4: per process work using the above grid data split
     pool = multiprocessing.Pool(num_cores)
     for sub_dictionary_of_cube_data in list_dictionaries_for_cores:
-        pool.apply_async(per_core_work, args=(
-            start_frame_coord_array, end_frame_coord_array, sub_dictionary_of_cube_data, MDA_selection, start_frame,
-            end_frame), callback=log_result_to_parent)
+        pool.apply_async(
+            per_core_work,
+            args=(
+                start_frame_coord_array,
+                end_frame_coord_array,
+                sub_dictionary_of_cube_data,
+                MDA_selection,
+                start_frame,
+                end_frame,
+            ),
+            callback=log_result_to_parent,
+        )
     pool.close()
     pool.join()
-    #so, at this stage the parent process now has a single dictionary with all the cube objects updated from all
+    # so, at this stage the parent process now has a single dictionary with all the cube objects updated from all
     # available cores
-    #the 3D streamplot (i.e, mayavi flow() function) will require separate 3D np arrays for dx,dy,dz
-    #the shape of each 3D array will unfortunately have to match the mgrid data structure (bit of a pain): (
+    # the 3D streamplot (i.e, mayavi flow() function) will require separate 3D np arrays for dx,dy,dz
+    # the shape of each 3D array will unfortunately have to match the mgrid data structure (bit of a pain): (
     # num_sheets - 1, num_sheets - 1, cubes_per_column)
     cubes_per_sheet = int(float(total_cubes) / float(num_sheets - 1))
-    #produce dummy zero arrays for dx,dy,dz of the appropriate shape:
+    # produce dummy zero arrays for dx,dy,dz of the appropriate shape:
     dx_array = np.zeros(delta_array_shape)
     dy_array = np.zeros(delta_array_shape)
     dz_array = np.zeros(delta_array_shape)
-    #now use the parent cube dictionary to correctly substitute in dx,dy,dz values
+    # now use the parent cube dictionary to correctly substitute in dx,dy,dz values
     current_sheet = 0  # which is also the current row
     y_index_current_sheet = 0  # sub row
     z_index_current_column = 0  # column
     total_cubes_current_sheet = 0
     for cube_number in range(0, total_cubes):
-        dx_array[current_sheet, y_index_current_sheet, z_index_current_column] = parent_cube_dictionary[cube_number][
-            'dx']
-        dy_array[current_sheet, y_index_current_sheet, z_index_current_column] = parent_cube_dictionary[cube_number][
-            'dy']
-        dz_array[current_sheet, y_index_current_sheet, z_index_current_column] = parent_cube_dictionary[cube_number][
-            'dz']
+        dx_array[
+            current_sheet, y_index_current_sheet, z_index_current_column
+        ] = parent_cube_dictionary[cube_number]["dx"]
+        dy_array[
+            current_sheet, y_index_current_sheet, z_index_current_column
+        ] = parent_cube_dictionary[cube_number]["dy"]
+        dz_array[
+            current_sheet, y_index_current_sheet, z_index_current_column
+        ] = parent_cube_dictionary[cube_number]["dz"]
         z_index_current_column += 1
         total_cubes_current_sheet += 1
         if z_index_current_column == delta_array_shape[2]:
             # done building current y-column so iterate y value and reset z
             z_index_current_column = 0
             y_index_current_sheet += 1
-            if y_index_current_sheet == delta_array_shape[1]:  # current sheet is complete
+            if (
+                y_index_current_sheet == delta_array_shape[1]
+            ):  # current sheet is complete
                 current_sheet += 1
                 y_index_current_sheet = 0  # restart for new sheet
                 z_index_current_column = 0
diff --git a/package/pyproject.toml b/package/pyproject.toml
index 297ed6f3040..9eabfc13553 100644
--- a/package/pyproject.toml
+++ b/package/pyproject.toml
@@ -129,6 +129,9 @@ include = '''
 (
 tables\.py
 | due\.py
+| setup\.py
+| visualization/.*\.py
 )
 '''
+extend-exclude = '__pycache__'
 required-version = '24'
diff --git a/package/setup.py b/package/setup.py
index 0fc4bef4e2a..412087cdaf2 100755
--- a/package/setup.py
+++ b/package/setup.py
@@ -60,7 +60,7 @@
 # NOTE: keep in sync with MDAnalysis.__version__ in version.py
 RELEASE = "2.8.0-dev0"
 
-is_release = 'dev' not in RELEASE
+is_release = "dev" not in RELEASE
 
 # Handle cython modules
 try:
@@ -68,18 +68,22 @@
     # minimum cython version now set to 0.28 to match pyproject.toml
     import Cython
     from Cython.Build import cythonize
+
     cython_found = True
 
     required_version = "0.28"
     if not Version(Cython.__version__) >= Version(required_version):
         # We don't necessarily die here. Maybe we already have
         #  the cythonized '.c' files.
-        print("Cython version {0} was found but won't be used: version {1} "
-              "or greater is required because it offers a handy "
-              "parallelization module".format(
-               Cython.__version__, required_version))
+        print(
+            "Cython version {0} was found but won't be used: version {1} "
+            "or greater is required because it offers a handy "
+            "parallelization module".format(
+                Cython.__version__, required_version
+            )
+        )
         cython_found = False
-    cython_linetrace = bool(os.environ.get('CYTHON_TRACE_NOGIL', False))
+    cython_linetrace = bool(os.environ.get("CYTHON_TRACE_NOGIL", False))
 except ImportError:
     cython_found = False
     if not is_release:
@@ -88,9 +92,10 @@
         sys.exit(1)
     cython_linetrace = False
 
+
 def abspath(file):
-    return os.path.join(os.path.dirname(os.path.abspath(__file__)),
-                        file)
+    return os.path.join(os.path.dirname(os.path.abspath(__file__)), file)
+
 
 class Config(object):
     """Config wrapper class to get build options
@@ -109,31 +114,31 @@ class Config(object):
 
     """
 
-    def __init__(self, fname='setup.cfg'):
+    def __init__(self, fname="setup.cfg"):
         fname = abspath(fname)
         if os.path.exists(fname):
             self.config = configparser.ConfigParser()
             self.config.read(fname)
 
     def get(self, option_name, default=None):
-        environ_name = 'MDA_' + option_name.upper()
+        environ_name = "MDA_" + option_name.upper()
         if environ_name in os.environ:
             val = os.environ[environ_name]
-            if val.upper() in ('1', 'TRUE'):
+            if val.upper() in ("1", "TRUE"):
                 return True
-            elif val.upper() in ('0', 'FALSE'):
+            elif val.upper() in ("0", "FALSE"):
                 return False
             return val
         try:
-            option = self.config.get('options', option_name)
+            option = self.config.get("options", option_name)
             return option
         except configparser.NoOptionError:
             return default
 
 
 class MDAExtension(Extension, object):
-    """Derived class to cleanly handle setup-time (numpy) dependencies.
-    """
+    """Derived class to cleanly handle setup-time (numpy) dependencies."""
+
     # The only setup-time numpy dependency comes when setting up its
     #  include dir.
     # The actual numpy import and call can be delayed until after pip
@@ -151,7 +156,7 @@ def include_dirs(self):
         if not self._mda_include_dirs:
             for item in self._mda_include_dir_args:
                 try:
-                    self._mda_include_dirs.append(item()) #The numpy callable
+                    self._mda_include_dirs.append(item())  # The numpy callable
                 except TypeError:
                     item = abspath(item)
                     self._mda_include_dirs.append((item))
@@ -174,9 +179,13 @@ def get_numpy_include():
         import numpy as np
     except ImportError:
         print('*** package "numpy" not found ***')
-        print('MDAnalysis requires a version of NumPy (>=1.21.0), even for setup.')
-        print('Please get it from http://numpy.scipy.org/ or install it through '
-              'your package manager.')
+        print(
+            "MDAnalysis requires a version of NumPy (>=1.21.0), even for setup."
+        )
+        print(
+            "Please get it from http://numpy.scipy.org/ or install it through "
+            "your package manager."
+        )
         sys.exit(-1)
     return np.get_include()
 
@@ -184,26 +193,27 @@ def get_numpy_include():
 def hasfunction(cc, funcname, include=None, extra_postargs=None):
     # From http://stackoverflow.com/questions/
     #            7018879/disabling-output-when-compiling-with-distutils
-    tmpdir = tempfile.mkdtemp(prefix='hasfunction-')
+    tmpdir = tempfile.mkdtemp(prefix="hasfunction-")
     devnull = oldstderr = None
     try:
         try:
-            fname = os.path.join(tmpdir, 'funcname.c')
-            with open(fname, 'w') as f:
+            fname = os.path.join(tmpdir, "funcname.c")
+            with open(fname, "w") as f:
                 if include is not None:
-                    f.write('#include {0!s}\n'.format(include))
-                f.write('int main(void) {\n')
-                f.write('    {0!s};\n'.format(funcname))
-                f.write('}\n')
+                    f.write("#include {0!s}\n".format(include))
+                f.write("int main(void) {\n")
+                f.write("    {0!s};\n".format(funcname))
+                f.write("}\n")
             # Redirect stderr to /dev/null to hide any error messages
             # from the compiler.
             # This will have to be changed if we ever have to check
             # for a function on Windows.
-            devnull = open('/dev/null', 'w')
+            devnull = open("/dev/null", "w")
             oldstderr = os.dup(sys.stderr.fileno())
             os.dup2(devnull.fileno(), sys.stderr.fileno())
-            objects = cc.compile([fname], output_dir=tmpdir,
-                                 extra_postargs=extra_postargs)
+            objects = cc.compile(
+                [fname], output_dir=tmpdir, extra_postargs=extra_postargs
+            )
             cc.link_executable(objects, os.path.join(tmpdir, "a.out"))
         except Exception:
             return False
@@ -221,11 +231,15 @@ def detect_openmp():
     print("Attempting to autodetect OpenMP support... ", end="")
     compiler = new_compiler()
     customize_compiler(compiler)
-    compiler.add_library('gomp')
-    include = '<omp.h>'
-    extra_postargs = ['-fopenmp']
-    hasopenmp = hasfunction(compiler, 'omp_get_num_threads()', include=include,
-                            extra_postargs=extra_postargs)
+    compiler.add_library("gomp")
+    include = "<omp.h>"
+    extra_postargs = ["-fopenmp"]
+    hasopenmp = hasfunction(
+        compiler,
+        "omp_get_num_threads()",
+        include=include,
+        extra_postargs=extra_postargs,
+    )
     if hasopenmp:
         print("Compiler supports OpenMP")
     else:
@@ -238,12 +252,12 @@ def using_clang():
     compiler = new_compiler()
     customize_compiler(compiler)
     compiler_ver = getoutput("{0} -v".format(compiler.compiler[0]))
-    if 'Spack GCC' in compiler_ver:
+    if "Spack GCC" in compiler_ver:
         # when gcc toolchain is built from source with spack
         # using clang, the 'clang' string may be present in
         # the compiler metadata, but it is not clang
         is_clang = False
-    elif 'clang' in compiler_ver:
+    elif "clang" in compiler_ver:
         # by default, Apple will typically alias gcc to
         # clang, with some mention of 'clang' in the
         # metadata
@@ -255,196 +269,252 @@ def using_clang():
 
 def extensions(config):
     # usually (except coming from release tarball) cython files must be generated
-    use_cython = config.get('use_cython', default=cython_found)
-    use_openmp = config.get('use_openmp', default=True)
-    annotate_cython = config.get('annotate_cython', default=False)
-
-    extra_compile_args = ['-std=c11', '-O3', '-funroll-loops',
-                          '-fsigned-zeros'] # see #2722
+    use_cython = config.get("use_cython", default=cython_found)
+    use_openmp = config.get("use_openmp", default=True)
+    annotate_cython = config.get("annotate_cython", default=False)
+
+    extra_compile_args = [
+        "-std=c11",
+        "-O3",
+        "-funroll-loops",
+        "-fsigned-zeros",
+    ]  # see #2722
     define_macros = []
-    if config.get('debug_cflags', default=False):
-        extra_compile_args.extend(['-Wall', '-pedantic'])
-        define_macros.extend([('DEBUG', '1')])
+    if config.get("debug_cflags", default=False):
+        extra_compile_args.extend(["-Wall", "-pedantic"])
+        define_macros.extend([("DEBUG", "1")])
 
     # encore is sensitive to floating point accuracy, especially on non-x86
     # to avoid reducing optimisations on everything, we make a set of compile
     # args specific to encore see #2997 for an example of this.
-    encore_compile_args = [a for a in extra_compile_args if 'O3' not in a]
-    if platform.machine() == 'aarch64' or platform.machine() == 'ppc64le':
-        encore_compile_args.append('-O1')
+    encore_compile_args = [a for a in extra_compile_args if "O3" not in a]
+    if platform.machine() == "aarch64" or platform.machine() == "ppc64le":
+        encore_compile_args.append("-O1")
     else:
-        encore_compile_args.append('-O3')
+        encore_compile_args.append("-O3")
 
     # allow using custom c/c++ flags and architecture specific instructions.
     # This allows people to build optimized versions of MDAnalysis.
     # Do here so not included in encore
-    extra_cflags = config.get('extra_cflags', default=False)
+    extra_cflags = config.get("extra_cflags", default=False)
     if extra_cflags:
         flags = extra_cflags.split()
         extra_compile_args.extend(flags)
 
-    cpp_extra_compile_args = [a for a in extra_compile_args if 'std' not in a]
-    cpp_extra_compile_args.append('-std=c++11')
-    cpp_extra_link_args=[]
+    cpp_extra_compile_args = [a for a in extra_compile_args if "std" not in a]
+    cpp_extra_compile_args.append("-std=c++11")
+    cpp_extra_link_args = []
     # needed to specify c++ runtime library on OSX
-    if platform.system() == 'Darwin' and using_clang():
-        cpp_extra_compile_args.append('-stdlib=libc++')
-        cpp_extra_compile_args.append('-mmacosx-version-min=10.9')
-        cpp_extra_link_args.append('-stdlib=libc++')
-        cpp_extra_link_args.append('-mmacosx-version-min=10.9')
+    if platform.system() == "Darwin" and using_clang():
+        cpp_extra_compile_args.append("-stdlib=libc++")
+        cpp_extra_compile_args.append("-mmacosx-version-min=10.9")
+        cpp_extra_link_args.append("-stdlib=libc++")
+        cpp_extra_link_args.append("-mmacosx-version-min=10.9")
 
     # Needed for large-file seeking under 32bit systems (for xtc/trr indexing
     # and access).
     largefile_macros = [
-        ('_LARGEFILE_SOURCE', None),
-        ('_LARGEFILE64_SOURCE', None),
-        ('_FILE_OFFSET_BITS', '64')
+        ("_LARGEFILE_SOURCE", None),
+        ("_LARGEFILE64_SOURCE", None),
+        ("_FILE_OFFSET_BITS", "64"),
     ]
 
     has_openmp = detect_openmp()
 
     if use_openmp and not has_openmp:
-        print('No openmp compatible compiler found default to serial build.')
+        print("No openmp compatible compiler found default to serial build.")
 
-    parallel_args = ['-fopenmp'] if has_openmp and use_openmp else []
-    parallel_libraries = ['gomp'] if has_openmp and use_openmp else []
-    parallel_macros = [('PARALLEL', None)] if has_openmp and use_openmp else []
+    parallel_args = ["-fopenmp"] if has_openmp and use_openmp else []
+    parallel_libraries = ["gomp"] if has_openmp and use_openmp else []
+    parallel_macros = [("PARALLEL", None)] if has_openmp and use_openmp else []
 
     if use_cython:
-        print('Will attempt to use Cython.')
+        print("Will attempt to use Cython.")
         if not cython_found:
-            print("Couldn't find a Cython installation. "
-                  "Not recompiling cython extensions.")
+            print(
+                "Couldn't find a Cython installation. "
+                "Not recompiling cython extensions."
+            )
             use_cython = False
     else:
-        print('Will not attempt to use Cython.')
+        print("Will not attempt to use Cython.")
 
-    source_suffix = '.pyx' if use_cython else '.c'
-    cpp_source_suffix = '.pyx' if use_cython else '.cpp'
+    source_suffix = ".pyx" if use_cython else ".c"
+    cpp_source_suffix = ".pyx" if use_cython else ".cpp"
 
     # The callable is passed so that it is only evaluated at install time.
 
     include_dirs = [get_numpy_include]
     # Windows automatically handles math library linking
     # and will not build MDAnalysis if we try to specify one
-    if os.name == 'nt':
+    if os.name == "nt":
         mathlib = []
     else:
-        mathlib = ['m']
+        mathlib = ["m"]
 
     if cython_linetrace:
         extra_compile_args.append("-DCYTHON_TRACE_NOGIL")
         cpp_extra_compile_args.append("-DCYTHON_TRACE_NOGIL")
 
-    libdcd = MDAExtension('MDAnalysis.lib.formats.libdcd',
-                          ['MDAnalysis/lib/formats/libdcd' + source_suffix],
-                          include_dirs=include_dirs + ['MDAnalysis/lib/formats/include'],
-                          define_macros=define_macros,
-                          extra_compile_args=extra_compile_args)
-    distances = MDAExtension('MDAnalysis.lib.c_distances',
-                             ['MDAnalysis/lib/c_distances' + source_suffix],
-                             include_dirs=include_dirs + ['MDAnalysis/lib/include'],
-                             libraries=mathlib,
-                             define_macros=define_macros,
-                             extra_compile_args=extra_compile_args)
-    distances_omp = MDAExtension('MDAnalysis.lib.c_distances_openmp',
-                                 ['MDAnalysis/lib/c_distances_openmp' + source_suffix],
-                                 include_dirs=include_dirs + ['MDAnalysis/lib/include'],
-                                 libraries=mathlib + parallel_libraries,
-                                 define_macros=define_macros + parallel_macros,
-                                 extra_compile_args=parallel_args + extra_compile_args,
-                                 extra_link_args=parallel_args)
-    qcprot = MDAExtension('MDAnalysis.lib.qcprot',
-                          ['MDAnalysis/lib/qcprot' + source_suffix],
-                          include_dirs=include_dirs,
-                          define_macros=define_macros,
-                          extra_compile_args=extra_compile_args)
-    transformation = MDAExtension('MDAnalysis.lib._transformations',
-                                  ['MDAnalysis/lib/src/transformations/transformations.c'],
-                                  libraries=mathlib,
-                                  define_macros=define_macros,
-                                  include_dirs=include_dirs,
-                                  extra_compile_args=extra_compile_args)
-    libmdaxdr = MDAExtension('MDAnalysis.lib.formats.libmdaxdr',
-                             sources=['MDAnalysis/lib/formats/libmdaxdr' + source_suffix,
-                                      'MDAnalysis/lib/formats/src/xdrfile.c',
-                                      'MDAnalysis/lib/formats/src/xdrfile_xtc.c',
-                                      'MDAnalysis/lib/formats/src/xdrfile_trr.c',
-                                      'MDAnalysis/lib/formats/src/trr_seek.c',
-                                      'MDAnalysis/lib/formats/src/xtc_seek.c',
-                             ],
-                             include_dirs=include_dirs + ['MDAnalysis/lib/formats/include',
-                                                          'MDAnalysis/lib/formats'],
-                             define_macros=largefile_macros + define_macros,
-                             extra_compile_args=extra_compile_args)
-    util = MDAExtension('MDAnalysis.lib.formats.cython_util',
-                        sources=['MDAnalysis/lib/formats/cython_util' + source_suffix],
-                        include_dirs=include_dirs,
-                        define_macros=define_macros,
-                        extra_compile_args=extra_compile_args)
-    cutil = MDAExtension('MDAnalysis.lib._cutil',
-                         sources=['MDAnalysis/lib/_cutil' + cpp_source_suffix],
-                         language='c++',
-                         libraries=mathlib,
-                         include_dirs=include_dirs + ['MDAnalysis/lib/include'],
-                         define_macros=define_macros,
-                         extra_compile_args=cpp_extra_compile_args,
-                         extra_link_args= cpp_extra_link_args)
-    augment = MDAExtension('MDAnalysis.lib._augment',
-                         sources=['MDAnalysis/lib/_augment' + cpp_source_suffix],
-                         language='c++',
-                         include_dirs=include_dirs,
-                         define_macros=define_macros,
-                         extra_compile_args=cpp_extra_compile_args,
-                         extra_link_args= cpp_extra_link_args)
-    timestep = MDAExtension('MDAnalysis.coordinates.timestep',
-                         sources=['MDAnalysis/coordinates/timestep' + cpp_source_suffix],
-                         language='c++',
-                         include_dirs=include_dirs,
-                         define_macros=define_macros,
-                         extra_compile_args=cpp_extra_compile_args,
-                         extra_link_args= cpp_extra_link_args)
-
-
-    encore_utils = MDAExtension('MDAnalysis.analysis.encore.cutils',
-                                sources=['MDAnalysis/analysis/encore/cutils' + source_suffix],
-                                include_dirs=include_dirs,
-                                define_macros=define_macros,
-                                extra_compile_args=encore_compile_args)
-    ap_clustering = MDAExtension('MDAnalysis.analysis.encore.clustering.affinityprop',
-                                 sources=['MDAnalysis/analysis/encore/clustering/affinityprop' + source_suffix,
-                                          'MDAnalysis/analysis/encore/clustering/src/ap.c'],
-                                 include_dirs=include_dirs+['MDAnalysis/analysis/encore/clustering/include'],
-                                 libraries=mathlib,
-                                 define_macros=define_macros,
-                                 extra_compile_args=encore_compile_args)
-    spe_dimred = MDAExtension('MDAnalysis.analysis.encore.dimensionality_reduction.stochasticproxembed',
-                              sources=['MDAnalysis/analysis/encore/dimensionality_reduction/stochasticproxembed' + source_suffix,
-                                       'MDAnalysis/analysis/encore/dimensionality_reduction/src/spe.c'],
-                              include_dirs=include_dirs+['MDAnalysis/analysis/encore/dimensionality_reduction/include'],
-                              libraries=mathlib,
-                              define_macros=define_macros,
-                              extra_compile_args=encore_compile_args)
-    nsgrid = MDAExtension('MDAnalysis.lib.nsgrid',
-                             ['MDAnalysis/lib/nsgrid' + cpp_source_suffix],
-                             include_dirs=include_dirs + ['MDAnalysis/lib/include'],
-                             language='c++',
-                             define_macros=define_macros,
-                             extra_compile_args=cpp_extra_compile_args,
-                             extra_link_args= cpp_extra_link_args)
-    pre_exts = [libdcd, distances, distances_omp, qcprot,
-                transformation, libmdaxdr, util, encore_utils,
-                ap_clustering, spe_dimred, cutil, augment, nsgrid, timestep]
+    libdcd = MDAExtension(
+        "MDAnalysis.lib.formats.libdcd",
+        ["MDAnalysis/lib/formats/libdcd" + source_suffix],
+        include_dirs=include_dirs + ["MDAnalysis/lib/formats/include"],
+        define_macros=define_macros,
+        extra_compile_args=extra_compile_args,
+    )
+    distances = MDAExtension(
+        "MDAnalysis.lib.c_distances",
+        ["MDAnalysis/lib/c_distances" + source_suffix],
+        include_dirs=include_dirs + ["MDAnalysis/lib/include"],
+        libraries=mathlib,
+        define_macros=define_macros,
+        extra_compile_args=extra_compile_args,
+    )
+    distances_omp = MDAExtension(
+        "MDAnalysis.lib.c_distances_openmp",
+        ["MDAnalysis/lib/c_distances_openmp" + source_suffix],
+        include_dirs=include_dirs + ["MDAnalysis/lib/include"],
+        libraries=mathlib + parallel_libraries,
+        define_macros=define_macros + parallel_macros,
+        extra_compile_args=parallel_args + extra_compile_args,
+        extra_link_args=parallel_args,
+    )
+    qcprot = MDAExtension(
+        "MDAnalysis.lib.qcprot",
+        ["MDAnalysis/lib/qcprot" + source_suffix],
+        include_dirs=include_dirs,
+        define_macros=define_macros,
+        extra_compile_args=extra_compile_args,
+    )
+    transformation = MDAExtension(
+        "MDAnalysis.lib._transformations",
+        ["MDAnalysis/lib/src/transformations/transformations.c"],
+        libraries=mathlib,
+        define_macros=define_macros,
+        include_dirs=include_dirs,
+        extra_compile_args=extra_compile_args,
+    )
+    libmdaxdr = MDAExtension(
+        "MDAnalysis.lib.formats.libmdaxdr",
+        sources=[
+            "MDAnalysis/lib/formats/libmdaxdr" + source_suffix,
+            "MDAnalysis/lib/formats/src/xdrfile.c",
+            "MDAnalysis/lib/formats/src/xdrfile_xtc.c",
+            "MDAnalysis/lib/formats/src/xdrfile_trr.c",
+            "MDAnalysis/lib/formats/src/trr_seek.c",
+            "MDAnalysis/lib/formats/src/xtc_seek.c",
+        ],
+        include_dirs=include_dirs
+        + ["MDAnalysis/lib/formats/include", "MDAnalysis/lib/formats"],
+        define_macros=largefile_macros + define_macros,
+        extra_compile_args=extra_compile_args,
+    )
+    util = MDAExtension(
+        "MDAnalysis.lib.formats.cython_util",
+        sources=["MDAnalysis/lib/formats/cython_util" + source_suffix],
+        include_dirs=include_dirs,
+        define_macros=define_macros,
+        extra_compile_args=extra_compile_args,
+    )
+    cutil = MDAExtension(
+        "MDAnalysis.lib._cutil",
+        sources=["MDAnalysis/lib/_cutil" + cpp_source_suffix],
+        language="c++",
+        libraries=mathlib,
+        include_dirs=include_dirs + ["MDAnalysis/lib/include"],
+        define_macros=define_macros,
+        extra_compile_args=cpp_extra_compile_args,
+        extra_link_args=cpp_extra_link_args,
+    )
+    augment = MDAExtension(
+        "MDAnalysis.lib._augment",
+        sources=["MDAnalysis/lib/_augment" + cpp_source_suffix],
+        language="c++",
+        include_dirs=include_dirs,
+        define_macros=define_macros,
+        extra_compile_args=cpp_extra_compile_args,
+        extra_link_args=cpp_extra_link_args,
+    )
+    timestep = MDAExtension(
+        "MDAnalysis.coordinates.timestep",
+        sources=["MDAnalysis/coordinates/timestep" + cpp_source_suffix],
+        language="c++",
+        include_dirs=include_dirs,
+        define_macros=define_macros,
+        extra_compile_args=cpp_extra_compile_args,
+        extra_link_args=cpp_extra_link_args,
+    )
 
+    encore_utils = MDAExtension(
+        "MDAnalysis.analysis.encore.cutils",
+        sources=["MDAnalysis/analysis/encore/cutils" + source_suffix],
+        include_dirs=include_dirs,
+        define_macros=define_macros,
+        extra_compile_args=encore_compile_args,
+    )
+    ap_clustering = MDAExtension(
+        "MDAnalysis.analysis.encore.clustering.affinityprop",
+        sources=[
+            "MDAnalysis/analysis/encore/clustering/affinityprop"
+            + source_suffix,
+            "MDAnalysis/analysis/encore/clustering/src/ap.c",
+        ],
+        include_dirs=include_dirs
+        + ["MDAnalysis/analysis/encore/clustering/include"],
+        libraries=mathlib,
+        define_macros=define_macros,
+        extra_compile_args=encore_compile_args,
+    )
+    spe_dimred = MDAExtension(
+        "MDAnalysis.analysis.encore.dimensionality_reduction.stochasticproxembed",
+        sources=[
+            "MDAnalysis/analysis/encore/dimensionality_reduction/stochasticproxembed"
+            + source_suffix,
+            "MDAnalysis/analysis/encore/dimensionality_reduction/src/spe.c",
+        ],
+        include_dirs=include_dirs
+        + ["MDAnalysis/analysis/encore/dimensionality_reduction/include"],
+        libraries=mathlib,
+        define_macros=define_macros,
+        extra_compile_args=encore_compile_args,
+    )
+    nsgrid = MDAExtension(
+        "MDAnalysis.lib.nsgrid",
+        ["MDAnalysis/lib/nsgrid" + cpp_source_suffix],
+        include_dirs=include_dirs + ["MDAnalysis/lib/include"],
+        language="c++",
+        define_macros=define_macros,
+        extra_compile_args=cpp_extra_compile_args,
+        extra_link_args=cpp_extra_link_args,
+    )
+    pre_exts = [
+        libdcd,
+        distances,
+        distances_omp,
+        qcprot,
+        transformation,
+        libmdaxdr,
+        util,
+        encore_utils,
+        ap_clustering,
+        spe_dimred,
+        cutil,
+        augment,
+        nsgrid,
+        timestep,
+    ]
 
     cython_generated = []
     if use_cython:
         extensions = cythonize(
             pre_exts,
             annotate=annotate_cython,
-            compiler_directives={'linetrace': cython_linetrace,
-                                 'embedsignature': False,
-                                 'language_level': '3'},
+            compiler_directives={
+                "linetrace": cython_linetrace,
+                "embedsignature": False,
+                "language_level": "3",
+            },
         )
         if cython_linetrace:
             print("Cython coverage will be enabled")
@@ -453,15 +523,16 @@ def extensions(config):
                 if source not in pre_ext.sources:
                     cython_generated.append(source)
     else:
-        #Let's check early for missing .c files
+        # Let's check early for missing .c files
         extensions = pre_exts
         for ext in extensions:
             for source in ext.sources:
-                if not (os.path.isfile(source) and
-                        os.access(source, os.R_OK)):
-                    raise IOError("Source file '{}' not found. This might be "
-                                "caused by a missing Cython install, or a "
-                                "failed/disabled Cython build.".format(source))
+                if not (os.path.isfile(source) and os.access(source, os.R_OK)):
+                    raise IOError(
+                        "Source file '{}' not found. This might be "
+                        "caused by a missing Cython install, or a "
+                        "failed/disabled Cython build.".format(source)
+                    )
     return extensions, cython_generated
 
 
@@ -477,7 +548,7 @@ def dynamic_author_list():
     "Chronological list of authors" title.
     """
     authors = []
-    with codecs.open(abspath('AUTHORS'), encoding='utf-8') as infile:
+    with codecs.open(abspath("AUTHORS"), encoding="utf-8") as infile:
         # An author is a bullet point under the title "Chronological list of
         # authors". We first want move the cursor down to the title of
         # interest.
@@ -486,21 +557,23 @@ def dynamic_author_list():
                 break
         else:
             # If we did not break, it means we did not find the authors.
-            raise IOError('EOF before the list of authors')
+            raise IOError("EOF before the list of authors")
         # Skip the next line as it is the title underlining
         line = next(infile)
         line_no += 1
-        if line[:4] != '----':
-            raise IOError('Unexpected content on line {0}, '
-                          'should be a string of "-".'.format(line_no))
+        if line[:4] != "----":
+            raise IOError(
+                "Unexpected content on line {0}, "
+                'should be a string of "-".'.format(line_no)
+            )
         # Add each bullet point as an author until the next title underlining
         for line in infile:
-            if line[:4] in ('----', '====', '~~~~'):
+            if line[:4] in ("----", "====", "~~~~"):
                 # The previous line was a title, hopefully it did not start as
                 # a bullet point so it got ignored. Since we hit a title, we
                 # are done reading the list of authors.
                 break
-            elif line.strip()[:2] == '- ':
+            elif line.strip()[:2] == "- ":
                 # This is a bullet point, so it should be an author name.
                 name = line.strip()[2:].strip()
                 authors.append(name)
@@ -509,28 +582,32 @@ def dynamic_author_list():
     # sorted alphabetically of the last name.
     authors.sort(key=lambda name: name.split()[-1])
     # Move Naveen and Elizabeth first, and Oliver last.
-    authors.remove('Naveen Michaud-Agrawal')
-    authors.remove('Elizabeth J. Denning')
-    authors.remove('Oliver Beckstein')
-    authors = (['Naveen Michaud-Agrawal', 'Elizabeth J. Denning']
-               + authors + ['Oliver Beckstein'])
+    authors.remove("Naveen Michaud-Agrawal")
+    authors.remove("Elizabeth J. Denning")
+    authors.remove("Oliver Beckstein")
+    authors = (
+        ["Naveen Michaud-Agrawal", "Elizabeth J. Denning"]
+        + authors
+        + ["Oliver Beckstein"]
+    )
 
     # Write the authors.py file.
-    out_path = abspath('MDAnalysis/authors.py')
-    with codecs.open(out_path, 'w', encoding='utf-8') as outfile:
+    out_path = abspath("MDAnalysis/authors.py")
+    with codecs.open(out_path, "w", encoding="utf-8") as outfile:
         # Write the header
-        header = '''\
+        header = """\
 #-*- coding:utf-8 -*-
 
 # This file is generated from the AUTHORS file during the installation process.
 # Do not edit it as your changes will be overwritten.
-'''
+"""
         print(header, file=outfile)
 
         # Write the list of authors as a python list
-        template = u'__authors__ = [\n{}\n]'
-        author_string = u',\n'.join(u'    u"{}"'.format(name)
-                                    for name in authors)
+        template = "__authors__ = [\n{}\n]"
+        author_string = ",\n".join(
+            '    u"{}"'.format(name) for name in authors
+        )
         print(template.format(author_string), file=outfile)
 
 
@@ -540,17 +617,18 @@ def long_description(readme):
     with open(abspath(readme)) as summary:
         buffer = summary.read()
     # remove top heading that messes up pypi display
-    m = re.search('====*\n[^\n]*README[^\n]*\n=====*\n', buffer,
-                  flags=re.DOTALL)
+    m = re.search(
+        "====*\n[^\n]*README[^\n]*\n=====*\n", buffer, flags=re.DOTALL
+    )
     assert m, "README.rst does not contain a level-1 heading"
-    return buffer[m.end():]
+    return buffer[m.end() :]
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     try:
         dynamic_author_list()
     except (OSError, IOError):
-        warnings.warn('Cannot write the list of authors.')
+        warnings.warn("Cannot write the list of authors.")
 
     try:
         # when building from repository for creating the distribution
@@ -563,24 +641,30 @@ def long_description(readme):
     config = Config()
     exts, cythonfiles = extensions(config)
 
-    setup(name='MDAnalysis',
-          version=RELEASE,
-          long_description=LONG_DESCRIPTION,
-          long_description_content_type='text/x-rst',
-          # currently unused & may become obsolte see setuptools #1569
-          provides=['MDAnalysis'],
-          ext_modules=exts,
-          test_suite="MDAnalysisTests",
-          tests_require=[
-              'MDAnalysisTests=={0!s}'.format(RELEASE),  # same as this release!
-          ],
+    setup(
+        name="MDAnalysis",
+        version=RELEASE,
+        long_description=LONG_DESCRIPTION,
+        long_description_content_type="text/x-rst",
+        # currently unused & may become obsolte see setuptools #1569
+        provides=["MDAnalysis"],
+        ext_modules=exts,
+        test_suite="MDAnalysisTests",
+        tests_require=[
+            "MDAnalysisTests=={0!s}".format(RELEASE),  # same as this release!
+        ],
     )
 
     # Releases keep their cythonized stuff for shipping.
-    if not config.get('keep_cythonized', default=is_release) and not cython_linetrace:
+    if (
+        not config.get("keep_cythonized", default=is_release)
+        and not cython_linetrace
+    ):
         for cythonized in cythonfiles:
             try:
                 os.unlink(cythonized)
             except OSError as err:
-                print("Warning: failed to delete cythonized file {0}: {1}. "
-                    "Moving on.".format(cythonized, err.strerror))
+                print(
+                    "Warning: failed to delete cythonized file {0}: {1}. "
+                    "Moving on.".format(cythonized, err.strerror)
+                )