diff --git a/poetry.lock b/poetry.lock
index beb12597e8dee2dd14b720cd5350b784d486e5e0..dffbd83e785a00dfc3b603e2fa77c5127fdd3ad1 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -64,6 +64,30 @@ python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*"
 [package.dependencies]
 pytz = ">=2015.7"
 
+[[package]]
+name = "black"
+version = "22.1.0"
+description = "The uncompromising code formatter."
+category = "dev"
+optional = false
+python-versions = ">=3.6.2"
+
+[package.dependencies]
+click = ">=8.0.0"
+dataclasses = {version = ">=0.6", markers = "python_version < \"3.7\""}
+mypy-extensions = ">=0.4.3"
+pathspec = ">=0.9.0"
+platformdirs = ">=2"
+tomli = ">=1.1.0"
+typed-ast = {version = ">=1.4.2", markers = "python_version < \"3.8\" and implementation_name == \"cpython\""}
+typing-extensions = {version = ">=3.10.0.0", markers = "python_version < \"3.10\""}
+
+[package.extras]
+colorama = ["colorama (>=0.4.3)"]
+d = ["aiohttp (>=3.7.4)"]
+jupyter = ["ipython (>=7.8.0)", "tokenize-rt (>=3.2.0)"]
+uvloop = ["uvloop (>=0.15.2)"]
+
 [[package]]
 name = "certifi"
 version = "2021.10.8"
@@ -83,6 +107,18 @@ python-versions = ">=3.5.0"
 [package.extras]
 unicode_backport = ["unicodedata2"]
 
+[[package]]
+name = "click"
+version = "8.0.4"
+description = "Composable command line interface toolkit"
+category = "dev"
+optional = false
+python-versions = ">=3.6"
+
+[package.dependencies]
+colorama = {version = "*", markers = "platform_system == \"Windows\""}
+importlib-metadata = {version = "*", markers = "python_version < \"3.8\""}
+
 [[package]]
 name = "colorama"
 version = "0.4.4"
@@ -118,7 +154,7 @@ name = "cowsay"
 version = "4.0"
 description = "The famous cowsay for GNU/Linux is now available for python"
 category = "main"
-optional = false
+optional = true
 python-versions = "*"
 
 [[package]]
@@ -365,6 +401,14 @@ numpy = ">=1.6"
 pyparsing = "*"
 scipy = "*"
 
+[[package]]
+name = "mypy-extensions"
+version = "0.4.3"
+description = "Experimental type system extensions for programs checked with the mypy typechecker."
+category = "dev"
+optional = false
+python-versions = "*"
+
 [[package]]
 name = "myst-parser"
 version = "0.13.7"
@@ -406,6 +450,14 @@ python-versions = ">=3.6"
 [package.dependencies]
 pyparsing = ">=2.0.2,<3.0.5 || >3.0.5"
 
+[[package]]
+name = "pathspec"
+version = "0.9.0"
+description = "Utility library for gitignore style pattern matching of file paths."
+category = "dev"
+optional = false
+python-versions = "!=3.0.*,!=3.1.*,!=3.2.*,!=3.3.*,!=3.4.*,>=2.7"
+
 [[package]]
 name = "pep8-naming"
 version = "0.10.0"
@@ -924,6 +976,14 @@ category = "dev"
 optional = false
 python-versions = ">=2.6, !=3.0.*, !=3.1.*, !=3.2.*"
 
+[[package]]
+name = "tomli"
+version = "1.2.3"
+description = "A lil' TOML parser"
+category = "dev"
+optional = false
+python-versions = ">=3.6"
+
 [[package]]
 name = "tox"
 version = "3.24.5"
@@ -1023,14 +1083,6 @@ category = "main"
 optional = false
 python-versions = "!=3.0.*,!=3.1.*,!=3.2.*,!=3.3.*,!=3.4.*,>=2.7"
 
-[[package]]
-name = "yapf"
-version = "0.31.0"
-description = "A formatter for Python code."
-category = "dev"
-optional = false
-python-versions = "*"
-
 [[package]]
 name = "zipp"
 version = "3.6.0"
@@ -1049,8 +1101,8 @@ docs = ["Sphinx", "sphinx-autoapi", "sphinx-rtd-theme", "myst-parser"]
 
 [metadata]
 lock-version = "1.1"
-python-versions = "^3.6"
-content-hash = "5f56853204a0856745eb9a253e385e0e78aad5aae07beeb448223306065f6b15"
+python-versions = "^3.6.2"
+content-hash = "bcf4d72f6cd263cc2d0668b472751140f6934655d2a34477d026089e6eff50bf"
 
 [metadata.files]
 alabaster = [
@@ -1077,6 +1129,31 @@ babel = [
     {file = "Babel-2.9.1-py2.py3-none-any.whl", hash = "sha256:ab49e12b91d937cd11f0b67cb259a57ab4ad2b59ac7a3b41d6c06c0ac5b0def9"},
     {file = "Babel-2.9.1.tar.gz", hash = "sha256:bc0c176f9f6a994582230df350aa6e05ba2ebe4b3ac317eab29d9be5d2768da0"},
 ]
+black = [
+    {file = "black-22.1.0-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:1297c63b9e1b96a3d0da2d85d11cd9bf8664251fd69ddac068b98dc4f34f73b6"},
+    {file = "black-22.1.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:2ff96450d3ad9ea499fc4c60e425a1439c2120cbbc1ab959ff20f7c76ec7e866"},
+    {file = "black-22.1.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:0e21e1f1efa65a50e3960edd068b6ae6d64ad6235bd8bfea116a03b21836af71"},
+    {file = "black-22.1.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:e2f69158a7d120fd641d1fa9a921d898e20d52e44a74a6fbbcc570a62a6bc8ab"},
+    {file = "black-22.1.0-cp310-cp310-win_amd64.whl", hash = "sha256:228b5ae2c8e3d6227e4bde5920d2fc66cc3400fde7bcc74f480cb07ef0b570d5"},
+    {file = "black-22.1.0-cp36-cp36m-macosx_10_9_x86_64.whl", hash = "sha256:b1a5ed73ab4c482208d20434f700d514f66ffe2840f63a6252ecc43a9bc77e8a"},
+    {file = "black-22.1.0-cp36-cp36m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:35944b7100af4a985abfcaa860b06af15590deb1f392f06c8683b4381e8eeaf0"},
+    {file = "black-22.1.0-cp36-cp36m-win_amd64.whl", hash = "sha256:7835fee5238fc0a0baf6c9268fb816b5f5cd9b8793423a75e8cd663c48d073ba"},
+    {file = "black-22.1.0-cp37-cp37m-macosx_10_9_x86_64.whl", hash = "sha256:dae63f2dbf82882fa3b2a3c49c32bffe144970a573cd68d247af6560fc493ae1"},
+    {file = "black-22.1.0-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:5fa1db02410b1924b6749c245ab38d30621564e658297484952f3d8a39fce7e8"},
+    {file = "black-22.1.0-cp37-cp37m-win_amd64.whl", hash = "sha256:c8226f50b8c34a14608b848dc23a46e5d08397d009446353dad45e04af0c8e28"},
+    {file = "black-22.1.0-cp38-cp38-macosx_10_9_universal2.whl", hash = "sha256:2d6f331c02f0f40aa51a22e479c8209d37fcd520c77721c034517d44eecf5912"},
+    {file = "black-22.1.0-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:742ce9af3086e5bd07e58c8feb09dbb2b047b7f566eb5f5bc63fd455814979f3"},
+    {file = "black-22.1.0-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:fdb8754b453fb15fad3f72cd9cad3e16776f0964d67cf30ebcbf10327a3777a3"},
+    {file = "black-22.1.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:f5660feab44c2e3cb24b2419b998846cbb01c23c7fe645fee45087efa3da2d61"},
+    {file = "black-22.1.0-cp38-cp38-win_amd64.whl", hash = "sha256:6f2f01381f91c1efb1451998bd65a129b3ed6f64f79663a55fe0e9b74a5f81fd"},
+    {file = "black-22.1.0-cp39-cp39-macosx_10_9_universal2.whl", hash = "sha256:efbadd9b52c060a8fc3b9658744091cb33c31f830b3f074422ed27bad2b18e8f"},
+    {file = "black-22.1.0-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:8871fcb4b447206904932b54b567923e5be802b9b19b744fdff092bd2f3118d0"},
+    {file = "black-22.1.0-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:ccad888050f5393f0d6029deea2a33e5ae371fd182a697313bdbd835d3edaf9c"},
+    {file = "black-22.1.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:07e5c049442d7ca1a2fc273c79d1aecbbf1bc858f62e8184abe1ad175c4f7cc2"},
+    {file = "black-22.1.0-cp39-cp39-win_amd64.whl", hash = "sha256:373922fc66676133ddc3e754e4509196a8c392fec3f5ca4486673e685a421321"},
+    {file = "black-22.1.0-py3-none-any.whl", hash = "sha256:3524739d76b6b3ed1132422bf9d82123cd1705086723bc3e235ca39fd21c667d"},
+    {file = "black-22.1.0.tar.gz", hash = "sha256:a7c0192d35635f6fc1174be575cb7915e92e5dd629ee79fdaf0dcfa41a80afb5"},
+]
 certifi = [
     {file = "certifi-2021.10.8-py2.py3-none-any.whl", hash = "sha256:d62a0163eb4c2344ac042ab2bdf75399a71a2d8c7d47eac2e2ee91b9d6339569"},
     {file = "certifi-2021.10.8.tar.gz", hash = "sha256:78884e7c1d4b00ce3cea67b44566851c4343c120abd683433ce934a68ea58872"},
@@ -1085,6 +1162,10 @@ charset-normalizer = [
     {file = "charset-normalizer-2.0.12.tar.gz", hash = "sha256:2857e29ff0d34db842cd7ca3230549d1a697f96ee6d3fb071cfa6c7393832597"},
     {file = "charset_normalizer-2.0.12-py3-none-any.whl", hash = "sha256:6881edbebdb17b39b4eaaa821b438bf6eddffb4468cf344f09f89def34a8b1df"},
 ]
+click = [
+    {file = "click-8.0.4-py3-none-any.whl", hash = "sha256:6a7a62563bbfabfda3a38f3023a1db4a35978c0abd76f6c9605ecd6554d6d9b1"},
+    {file = "click-8.0.4.tar.gz", hash = "sha256:8458d7b1287c5fb128c90e23381cf99dcde74beaf6c7ff6384ce84d6fe090adb"},
+]
 colorama = [
     {file = "colorama-0.4.4-py2.py3-none-any.whl", hash = "sha256:9f47eda37229f68eee03b24b9748937c7dc3868f906e8ba69fbcbdd3bc5dc3e2"},
     {file = "colorama-0.4.4.tar.gz", hash = "sha256:5941b2b48a20143d2267e95b1c2a7603ce057ee39fd88e7329b0c292aa16869b"},
@@ -1319,6 +1400,10 @@ mdplus = [
 mdtraj = [
     {file = "mdtraj-1.9.7.tar.gz", hash = "sha256:8a3309d2ef6ddb1023dcf48300d5df9b190469b63f69af9d55490bc4799d3757"},
 ]
+mypy-extensions = [
+    {file = "mypy_extensions-0.4.3-py2.py3-none-any.whl", hash = "sha256:090fedd75945a69ae91ce1303b5824f428daf5a028d2f6ab8a299250a846f15d"},
+    {file = "mypy_extensions-0.4.3.tar.gz", hash = "sha256:2d82818f5bb3e369420cb3c4060a7970edba416647068eb4c5343488a6c604a8"},
+]
 myst-parser = [
     {file = "myst-parser-0.13.7.tar.gz", hash = "sha256:e4bc99e43e19f70d22e528de8e7cce59f7e8e7c4c34dcba203de92de7a7c7c85"},
     {file = "myst_parser-0.13.7-py3-none-any.whl", hash = "sha256:260355b4da8e8865fe080b0638d7f1ab1791dc4bed02a7a48630b6bad4249219"},
@@ -1357,6 +1442,10 @@ packaging = [
     {file = "packaging-21.3-py3-none-any.whl", hash = "sha256:ef103e05f519cdc783ae24ea4e2e0f508a9c99b2d4969652eed6a2e1ea5bd522"},
     {file = "packaging-21.3.tar.gz", hash = "sha256:dd47c42927d89ab911e606518907cc2d3a1f38bbd026385970643f9c5b8ecfeb"},
 ]
+pathspec = [
+    {file = "pathspec-0.9.0-py2.py3-none-any.whl", hash = "sha256:7d15c4ddb0b5c802d161efc417ec1a2558ea2653c2e8ad9c19098201dc1c993a"},
+    {file = "pathspec-0.9.0.tar.gz", hash = "sha256:e564499435a2673d586f6b2130bb5b95f04a3ba06f81b8f895b651a3c76aabb1"},
+]
 pep8-naming = [
     {file = "pep8-naming-0.10.0.tar.gz", hash = "sha256:f3b4a5f9dd72b991bf7d8e2a341d2e1aa3a884a769b5aaac4f56825c1763bf3a"},
     {file = "pep8_naming-0.10.0-py2.py3-none-any.whl", hash = "sha256:5d9f1056cb9427ce344e98d1a7f5665710e2f20f748438e308995852cfa24164"},
@@ -1594,6 +1683,10 @@ toml = [
     {file = "toml-0.10.2-py2.py3-none-any.whl", hash = "sha256:806143ae5bfb6a3c6e736a764057db0e6a0e05e338b5630894a5f779cabb4f9b"},
     {file = "toml-0.10.2.tar.gz", hash = "sha256:b3bda1d108d5dd99f4a20d24d9c348e91c4db7ab1b749200bded2f839ccbe68f"},
 ]
+tomli = [
+    {file = "tomli-1.2.3-py3-none-any.whl", hash = "sha256:e3069e4be3ead9668e21cb9b074cd948f7b3113fd9c8bba083f48247aab8b11c"},
+    {file = "tomli-1.2.3.tar.gz", hash = "sha256:05b6166bff487dc068d322585c7ea4ef78deed501cc124060e0f238e89a9231f"},
+]
 tox = [
     {file = "tox-3.24.5-py2.py3-none-any.whl", hash = "sha256:be3362472a33094bce26727f5f771ca0facf6dafa217f65875314e9a6600c95c"},
     {file = "tox-3.24.5.tar.gz", hash = "sha256:67e0e32c90e278251fea45b696d0fef3879089ccbe979b0c556d35d5a70e2993"},
@@ -1698,10 +1791,6 @@ wrapt = [
     {file = "wrapt-1.13.3-cp39-cp39-win_amd64.whl", hash = "sha256:81bd7c90d28a4b2e1df135bfbd7c23aee3050078ca6441bead44c42483f9ebfb"},
     {file = "wrapt-1.13.3.tar.gz", hash = "sha256:1fea9cd438686e6682271d36f3481a9f3636195578bab9ca3382e2f5f01fc185"},
 ]
-yapf = [
-    {file = "yapf-0.31.0-py2.py3-none-any.whl", hash = "sha256:e3a234ba8455fe201eaa649cdac872d590089a18b661e39bbac7020978dd9c2e"},
-    {file = "yapf-0.31.0.tar.gz", hash = "sha256:408fb9a2b254c302f49db83c59f9aa0b4b0fd0ec25be3a5c51181327922ff63d"},
-]
 zipp = [
     {file = "zipp-3.6.0-py3-none-any.whl", hash = "sha256:9fe5ea21568a0a70e50f273397638d39b03353731e6cbbb3fd8502a33fec40bc"},
     {file = "zipp-3.6.0.tar.gz", hash = "sha256:71c644c5369f4a6e07636f0aa966270449561fcea2e3d6747b8d23efaa9d7832"},
diff --git a/pycgtool/__main__.py b/pycgtool/__main__.py
index ae46fed1a6b14793baa5af90de659ca2d987e25d..fde95790a7a808910b0dad14f58a698bbf54b458 100755
--- a/pycgtool/__main__.py
+++ b/pycgtool/__main__.py
@@ -34,7 +34,8 @@ class PyCGTOOL:
             topology_file=self.config.topology,
             trajectory_file=self.config.trajectory,  # May be None
             frame_start=self.config.begin,
-            frame_end=self.config.end)
+            frame_end=self.config.end,
+        )
 
         self.mapping = None
         self.out_frame = self.in_frame
@@ -43,7 +44,9 @@ class PyCGTOOL:
 
         if self.out_frame.natoms == 0:
             # The following steps won't work with frames containing no atoms
-            logger.warning('Mapped trajectory contains no atoms - check your mapping file is correct!')
+            logger.warning(
+                "Mapped trajectory contains no atoms - check your mapping file is correct!"
+            )
             return
 
         self.bondset = None
@@ -52,7 +55,7 @@ class PyCGTOOL:
             self.measure_bonds()
 
         if self.config.output_xtc:
-            self.out_frame.save(self.get_output_filepath('xtc'))
+            self.out_frame.save(self.get_output_filepath("xtc"))
 
     def get_output_filepath(self, ext: PathLike) -> pathlib.Path:
         """Get file path for an output file by extension.
@@ -61,22 +64,22 @@ class PyCGTOOL:
         """
         out_dir = pathlib.Path(self.config.out_dir)
         out_dir.mkdir(parents=True, exist_ok=True)
-        return out_dir.joinpath(self.config.output_name + '.' + ext)
+        return out_dir.joinpath(self.config.output_name + "." + ext)
 
     def apply_mapping(self, in_frame: Frame) -> typing.Tuple[Mapping, Frame]:
         """Map input frame to output using requested mapping file."""
-        mapping = Mapping(self.config.mapping,
-                          self.config,
-                          itp_filename=self.config.itp)
+        mapping = Mapping(
+            self.config.mapping, self.config, itp_filename=self.config.itp
+        )
         out_frame = mapping.apply(in_frame)
-        out_frame.save(self.get_output_filepath('gro'), frame_number=0)
+        out_frame.save(self.get_output_filepath("gro"), frame_number=0)
 
         if self.config.backmapper_resname and self.out_frame.n_frames > 1:
             try:
                 self.train_backmapper(self.config.resname)
 
             except ImportError:
-                logger.error('MDPlus must be installed to perform backmapping')
+                logger.error("MDPlus must be installed to perform backmapping")
 
         return mapping, out_frame
 
@@ -87,46 +90,47 @@ class PyCGTOOL:
         if self.mapping is not None and self.out_frame.n_frames > 1:
             # Only perform Boltzmann Inversion if we have a mapping and a trajectory.
             # Otherwise we get infinite force constants.
-            logger.info('Starting Boltzmann Inversion')
+            logger.info("Starting Boltzmann Inversion")
             self.bondset.boltzmann_invert()
-            logger.info('Finished Boltzmann Inversion')
+            logger.info("Finished Boltzmann Inversion")
 
             if self.config.output_forcefield:
                 logger.info("Writing GROMACS forcefield directory")
                 out_dir = pathlib.Path(self.config.out_dir)
-                forcefield = ForceField(self.config.output_name,
-                                        dir_path=out_dir)
-                forcefield.write(self.config.output_name, self.mapping,
-                                 self.bondset)
+                forcefield = ForceField(self.config.output_name, dir_path=out_dir)
+                forcefield.write(self.config.output_name, self.mapping, self.bondset)
                 logger.info("Finished writing GROMACS forcefield directory")
 
             else:
-                self.bondset.write_itp(self.get_output_filepath('itp'),
-                                       mapping=self.mapping)
+                self.bondset.write_itp(
+                    self.get_output_filepath("itp"), mapping=self.mapping
+                )
 
         if self.config.dump_measurements:
-            logger.info('Writing bond measurements to file')
-            self.bondset.dump_values(self.config.dump_n_values,
-                                     self.config.out_dir)
-            logger.info('Finished writing bond measurements to file')
+            logger.info("Writing bond measurements to file")
+            self.bondset.dump_values(self.config.dump_n_values, self.config.out_dir)
+            logger.info("Finished writing bond measurements to file")
 
     def train_backmapper(self, resname: str):
         from mdplus.multiscale import GLIMPS
-        sel = f'resname {resname}'
+
+        sel = f"resname {resname}"
 
         aa_subset_traj = self.in_frame._trajectory.atom_slice(
-            self.in_frame._trajectory.topology.select(sel))
+            self.in_frame._trajectory.topology.select(sel)
+        )
         cg_subset_traj = self.out_frame._trajectory.atom_slice(
-            self.out_frame._trajectory.topology.select(sel))
+            self.out_frame._trajectory.topology.select(sel)
+        )
 
-        logger.info('Training backmapper')
+        logger.info("Training backmapper")
         # Param x_valence is approximate number of bonds per CG bead
         # Values greater than 2 fail for small molecules e.g. sugar test case
         backmapper = GLIMPS(x_valence=2)
         backmapper.fit(cg_subset_traj.xyz, aa_subset_traj.xyz)
-        logger.info('Finished training backmapper')
+        logger.info("Finished training backmapper")
 
-        backmapper.save(self.get_output_filepath('backmapper.pkl'))
+        backmapper.save(self.get_output_filepath("backmapper.pkl"))
 
 
 class BooleanAction(argparse.Action):
@@ -134,15 +138,16 @@ class BooleanAction(argparse.Action):
 
     Based on https://thisdataguy.com/2017/07/03/no-options-with-argparse-and-python/
     """
+
     def __init__(self, option_strings, dest, nargs=None, **kwargs):
         super().__init__(option_strings, dest, nargs=0, **kwargs)
 
     def __call__(self, parser, namespace, values, option_string=None):
-        setattr(namespace, self.dest, not option_string.startswith('--no-'))
+        setattr(namespace, self.dest, not option_string.startswith("--no-"))
 
 
 def parse_arguments(arg_list):
-    # yapf: disable
+    # fmt: off
     parser = argparse.ArgumentParser(
         prog='pycgtool',
         formatter_class=argparse.ArgumentDefaultsHelpFormatter,
@@ -236,7 +241,7 @@ def parse_arguments(arg_list):
                                 help=argparse.SUPPRESS)
     secret_options.add_argument('--cow', default=False, action='store_true',
                                 help=argparse.SUPPRESS)
-    # yapf: enable
+    # fmt: on
 
     args = parser.parse_args(arg_list)
 
@@ -257,15 +262,15 @@ def validate_arguments(args):
     if not args.dump_measurements:
         args.dump_measurements = bool(args.bondset) and not bool(args.mapping)
         logger.info(
-            'Argument --dump-measurements has been set because you have provided a bondset but no mapping'
+            "Argument --dump-measurements has been set because you have provided a bondset but no mapping"
         )
 
     if not args.mapping and not args.bondset:
-        raise ArgumentValidationError('One or both of -m and -b is required')
+        raise ArgumentValidationError("One or both of -m and -b is required")
 
     if args.backmapper_resname:
         logger.warning(
-            'Backmapping is an experimental feature and has not yet been fully validated'
+            "Backmapping is an experimental feature and has not yet been fully validated"
         )
 
     return args
@@ -279,10 +284,12 @@ def main():
         logging.basicConfig(level=args.log_level)
 
     else:
-        logging.basicConfig(level=args.log_level,
-                            format='%(message)s',
-                            datefmt='[%X]',
-                            handlers=[RichHandler(rich_tracebacks=True)])
+        logging.basicConfig(
+            level=args.log_level,
+            format="%(message)s",
+            datefmt="[%X]",
+            handlers=[RichHandler(rich_tracebacks=True)],
+        )
 
         banner = r"""
              _____        _____ _____ _______ ____   ____  _      
@@ -299,37 +306,39 @@ def main():
         if args.cow:
             try:
                 import cowsay
-                banner = cowsay.cow('PyCGTOOL')
+
+                banner = cowsay.cow("PyCGTOOL")
 
             except ImportError:
                 pass
 
         else:
-            logger.info('[bold blue]%s[/]',
-                        textwrap.dedent(banner),
-                        extra={'markup': True})
+            logger.info(
+                "[bold blue]%s[/]", textwrap.dedent(banner), extra={"markup": True}
+            )
 
-    logger.info(30 * '-')
-    logger.info('Topology:\t%s', args.topology)
-    logger.info('Trajectory:\t%s', args.trajectory)
-    logger.info('Mapping:\t%s', args.mapping)
-    logger.info('Bondset:\t%s', args.bondset)
-    logger.info(30 * '-')
+    logger.info(30 * "-")
+    logger.info("Topology:\t%s", args.topology)
+    logger.info("Trajectory:\t%s", args.trajectory)
+    logger.info("Mapping:\t%s", args.mapping)
+    logger.info("Bondset:\t%s", args.bondset)
+    logger.info(30 * "-")
 
     try:
         if args.profile:
             with cProfile.Profile() as profiler:
                 pycgtool = PyCGTOOL(args)
 
-            profiler.dump_stats('gprof.out')
+            profiler.dump_stats("gprof.out")
 
         else:
             pycgtool = PyCGTOOL(args)
 
         elapsed_time = time.time() - start_time
-        logger.info('Processed %d frames in %.2f s',
-                    pycgtool.out_frame.n_frames, elapsed_time)
-        logger.info('Finished processing - goodbye!')
+        logger.info(
+            "Processed %d frames in %.2f s", pycgtool.out_frame.n_frames, elapsed_time
+        )
+        logger.info("Finished processing - goodbye!")
 
     except Exception as exc:
         logger.exception(exc)
diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py
index 4910b7d90291d431799c4a18b2134e7a5ad89b21..3a4dca5a41f8236b10ba3ae387567c5c85bb1efa 100644
--- a/pycgtool/bondset.py
+++ b/pycgtool/bondset.py
@@ -34,10 +34,13 @@ class Bond:
 
     Bond lengths, angles and dihedrals are all equivalent, distinguished by the number of atoms present.
     """
-    def __init__(self,
-                 atoms: typing.Iterable[str],
-                 atom_numbers: typing.Optional[typing.Iterable[int]] = None,
-                 func_form=None):
+
+    def __init__(
+        self,
+        atoms: typing.Iterable[str],
+        atom_numbers: typing.Optional[typing.Iterable[int]] = None,
+        func_form=None,
+    ):
         """Create a single bond definition.
 
         :param List[str] atoms: List of atom names defining the bond
@@ -62,8 +65,9 @@ class Bond:
         :param temp: Temperature at which the simulation was performed
         """
         if len(self.values) == 0:
-            raise ValueError("No bonds were measured between beads {0}".format(
-                self.atoms))
+            raise ValueError(
+                "No bonds were measured between beads {0}".format(self.atoms)
+            )
 
         values = self.values
         if not isinstance(self.values, np.ndarray):
@@ -77,7 +81,7 @@ class Bond:
             if len(self.atoms) == 4:
                 two_pi = 2 * np.pi
                 self.eqm -= np.pi
-                self.eqm = (((self.eqm + np.pi) % two_pi) - np.pi)
+                self.eqm = ((self.eqm + np.pi) % two_pi) - np.pi
 
             try:
                 self.fconst = self._func_form.fconst(values, temp)
@@ -89,7 +93,8 @@ class Bond:
     def __repr__(self) -> str:
         try:
             return "<Bond containing atoms {0} with r_0 {1:.3f} and force constant {2:.3e}>".format(
-                ", ".join(self.atoms), self.eqm, self.fconst)
+                ", ".join(self.atoms), self.eqm, self.fconst
+            )
         except (AttributeError, TypeError):
             return "<Bond containing atoms {0}>".format(", ".join(self.atoms))
 
@@ -99,6 +104,7 @@ class BondSet:
 
     BondSet contains a dictionary of lists of Bonds. Each list corresponds to a single molecule.
     """
+
     def _get_default_func_forms(self, options):
 
         try:
@@ -109,33 +115,32 @@ class BondSet:
 
         if default_fc:
             default_forms = [
-                'MartiniDefaultLength', 'MartiniDefaultAngle',
-                'MartiniDefaultDihedral'
+                "MartiniDefaultLength",
+                "MartiniDefaultAngle",
+                "MartiniDefaultDihedral",
             ]
         else:
-            default_forms = ['Harmonic', 'CosHarmonic', 'Harmonic']
+            default_forms = ["Harmonic", "CosHarmonic", "Harmonic"]
 
         functional_forms_map = get_functional_forms()
 
         functional_forms = [None, None]
         functional_forms.extend(
-            [functional_forms_map[ff].value for ff in default_forms])
+            [functional_forms_map[ff].value for ff in default_forms]
+        )
 
         try:
-            functional_forms[2] = functional_forms_map[
-                options.length_form].value
+            functional_forms[2] = functional_forms_map[options.length_form].value
         except AttributeError:
             pass
 
         try:
-            functional_forms[3] = functional_forms_map[
-                options.angle_form].value
+            functional_forms[3] = functional_forms_map[options.angle_form].value
         except AttributeError:
             pass
 
         try:
-            functional_forms[4] = functional_forms_map[
-                options.dihedral_form].value
+            functional_forms[4] = functional_forms_map[options.dihedral_form].value
         except AttributeError:
             pass
 
@@ -154,7 +159,7 @@ class BondSet:
         try:
             self._temperature = options.temperature
         except AttributeError:
-            self._temperature = 310.
+            self._temperature = 310.0
 
         self._functional_forms = self._get_default_func_forms(options)
 
@@ -166,19 +171,26 @@ class BondSet:
                 angles_defined = False
                 for atomlist in mol_section:
                     is_angle_or_dihedral = len(atomlist) > 2
-                    mean_function = circular_mean if is_angle_or_dihedral else np.nanmean
-                    variance_function = circular_variance if is_angle_or_dihedral else np.nanvar
+                    mean_function = (
+                        circular_mean if is_angle_or_dihedral else np.nanmean
+                    )
+                    variance_function = (
+                        circular_variance if is_angle_or_dihedral else np.nanvar
+                    )
 
                     # Construct instance of Boltzmann Inversion function and
                     # inject dependencies for mean and variance functions
                     # TODO: Should we allow overriding functional forms per bond?
                     func_form = self._functional_forms[len(atomlist)](
-                        mean_function, variance_function)
+                        mean_function, variance_function
+                    )
 
                     if {x for x in atomlist if atomlist.count(x) > 1}:
                         raise ValueError(
-                            "Defined bond '{0}' contains duplicate atoms".
-                            format(atomlist))
+                            "Defined bond '{0}' contains duplicate atoms".format(
+                                atomlist
+                            )
+                        )
 
                     mol_bonds.append(Bond(atoms=atomlist, func_form=func_form))
                     if len(atomlist) > 2:
@@ -190,16 +202,24 @@ class BondSet:
                     if options.generate_angles:
                         for atomlist in angles:
                             mol_bonds.append(
-                                Bond(atoms=atomlist,
-                                     func_form=self._functional_forms[3](
-                                         circular_mean, circular_variance)))
+                                Bond(
+                                    atoms=atomlist,
+                                    func_form=self._functional_forms[3](
+                                        circular_mean, circular_variance
+                                    ),
+                                )
+                            )
 
                     if options.generate_dihedrals:
                         for atomlist in dihedrals:
                             mol_bonds.append(
-                                Bond(atoms=atomlist,
-                                     func_form=self._functional_forms[4](
-                                         circular_mean, circular_variance)))
+                                Bond(
+                                    atoms=atomlist,
+                                    func_form=self._functional_forms[4](
+                                        circular_mean, circular_variance
+                                    ),
+                                )
+                            )
 
     @staticmethod
     def _create_angles(mol_bonds):
@@ -217,7 +237,7 @@ class BondSet:
         self,
         mol: str,
         natoms: int,
-        select: typing.Callable[[Bond], bool] = lambda x: True
+        select: typing.Callable[[Bond], bool] = lambda x: True,
     ) -> typing.List[Bond]:
         """Return list of bonds from molecule containing natoms atoms.
 
@@ -230,7 +250,8 @@ class BondSet:
             return [bond for bond in self._molecules[mol] if select(bond)]
 
         return [
-            bond for bond in self._molecules[mol]
+            bond
+            for bond in self._molecules[mol]
             if len(bond.atoms) == natoms and select(bond)
         ]
 
@@ -242,13 +263,12 @@ class BondSet:
         :return: List of bonds
         """
         if with_constr:
-            return [
-                bond for bond in self._molecules[mol] if len(bond.atoms) == 2
-            ]
+            return [bond for bond in self._molecules[mol] if len(bond.atoms) == 2]
 
         return [
-            bond for bond in self._molecules[mol] if len(bond.atoms) == 2
-            and bond.fconst < self._fconst_constr_threshold
+            bond
+            for bond in self._molecules[mol]
+            if len(bond.atoms) == 2 and bond.fconst < self._fconst_constr_threshold
         ]
 
     def get_bond_length_constraints(self, mol):
@@ -258,8 +278,9 @@ class BondSet:
         :return: List of bonds
         """
         return [
-            bond for bond in self._molecules[mol] if len(bond.atoms) == 2
-            and bond.fconst >= self._fconst_constr_threshold
+            bond
+            for bond in self._molecules[mol]
+            if len(bond.atoms) == 2 and bond.fconst >= self._fconst_constr_threshold
         ]
 
     def get_bond_angles(self, mol, exclude_triangle=True):
@@ -269,9 +290,7 @@ class BondSet:
         :param exclude_triangle: Exclude angles that are part of a triangle?
         :return: List of bonds
         """
-        angles = [
-            bond for bond in self._molecules[mol] if len(bond.atoms) == 3
-        ]
+        angles = [bond for bond in self._molecules[mol] if len(bond.atoms) == 3]
 
         if exclude_triangle:
             edges = [
@@ -282,15 +301,14 @@ class BondSet:
             def is_triangle(atoms):
                 triangle_edges = 0
                 for j in range(3):
-                    if (atoms[j - 1],
-                            atoms[j]) in edges or (atoms[j],
-                                                   atoms[j - 1]) in edges:
+                    if (atoms[j - 1], atoms[j]) in edges or (
+                        atoms[j],
+                        atoms[j - 1],
+                    ) in edges:
                         triangle_edges += 1
                 return triangle_edges >= 3
 
-            angles = [
-                angle for angle in angles if not is_triangle(angle.atoms)
-            ]
+            angles = [angle for angle in angles if not is_triangle(angle.atoms)]
 
         return angles
 
@@ -332,12 +350,13 @@ class BondSet:
                     ]
                 except ValueError as e:
                     missing = [
-                        atom for atom in bond.atoms
-                        if atom.lstrip("+-") not in index
+                        atom for atom in bond.atoms if atom.lstrip("+-") not in index
                     ]
                     e.args = (
                         "Bead(s) {0} do(es) not exist in residue {1}".format(
-                            missing, mol), )
+                            missing, mol
+                        ),
+                    )
                     raise
 
     def write_itp(self, filename, mapping):
@@ -352,22 +371,19 @@ class BondSet:
     def itp_text(self, mapping):
         atom_template = {
             "nomass": "{0:4d} {1:4s} {2:4d} {3:4s} {4:4s} {5:4d} {6:8.3f}",
-            "mass":
-            "{0:4d} {1:4s} {2:4d} {3:4s} {4:4s} {5:4d} {6:8.3f} {7:8.3f}"
+            "mass": "{0:4d} {1:4s} {2:4d} {3:4s} {4:4s} {5:4d} {6:8.3f} {7:8.3f}",
         }
 
-        def write_bond_angle_dih(bonds,
-                                 section_header,
-                                 print_fconst=True,
-                                 multiplicity=None,
-                                 rad2deg=False):
+        def write_bond_angle_dih(
+            bonds, section_header, print_fconst=True, multiplicity=None, rad2deg=False
+        ):
             ret_lines = []
             if bonds:
                 ret_lines.append("\n[ {0:s} ]".format(section_header))
             for bond in bonds:
-                line = " ".join([
-                    "{0:4d}".format(atnum + 1) for atnum in bond.atom_numbers
-                ])
+                line = " ".join(
+                    ["{0:4d}".format(atnum + 1) for atnum in bond.atom_numbers]
+                )
                 eqm = math.degrees(bond.eqm) if rad2deg else bond.eqm
                 line += " {0:4d} {1:12.5f}".format(bond.gromacs_type_id, eqm)
                 if print_fconst:
@@ -384,20 +400,29 @@ class BondSet:
             "; James Graham <J.A.Graham@soton.ac.uk> 2016",
             "; University of Southampton",
             "; https://github.com/jag1g13/pycgtool",
-            ";"
-        ]  # yapf: disable
+            ";",
+        ]  # fmt: skip
 
         # Print molecule
         for mol in self._molecules:
             if mol not in mapping:
-                logger.warning((
-                    "Molecule '%s' present in bonding file, but not in mapping. "
-                    "Parameters have not been calculated."), mol)
+                logger.warning(
+                    (
+                        "Molecule '%s' present in bonding file, but not in mapping. "
+                        "Parameters have not been calculated."
+                    ),
+                    mol,
+                )
                 continue
 
             if not all(bond.fconst is not None for bond in self._molecules[mol]):
-                logger.warning(("Molecule '%s' has no measured bond values. "
-                                "Parameters have not been calculated."), mol)
+                logger.warning(
+                    (
+                        "Molecule '%s' has no measured bond values. "
+                        "Parameters have not been calculated."
+                    ),
+                    mol,
+                )
                 continue
 
             ret_lines.append("\n[ moleculetype ]")
@@ -408,48 +433,63 @@ class BondSet:
             for i, bead in enumerate(mapping[mol], start=1):
                 # atnum type resnum resname atname c-group charge (mass)
                 if isinstance(bead, VirtualMap):
-                    ret_lines.append(atom_template["mass"].format(
-                        i, bead.type, 1, mol, bead.name, i, bead.charge,
-                        bead.mass))
+                    ret_lines.append(
+                        atom_template["mass"].format(
+                            i, bead.type, 1, mol, bead.name, i, bead.charge, bead.mass
+                        )
+                    )
                 else:
-                    ret_lines.append(atom_template["nomass"].format(
-                        i, bead.type, 1, mol, bead.name, i, bead.charge))
+                    ret_lines.append(
+                        atom_template["nomass"].format(
+                            i, bead.type, 1, mol, bead.name, i, bead.charge
+                        )
+                    )
 
             virtual_beads = self.get_virtual_beads(mapping[mol])
             if len(virtual_beads) != 0:
                 ret_lines.append("\n[ virtual_sitesn ]")
-                excl_lines = ["\n[ exclusions ]"
-                              ]  # Exclusions section for virtual sites
+                excl_lines = [
+                    "\n[ exclusions ]"
+                ]  # Exclusions section for virtual sites
 
                 for vbead in virtual_beads:
-                    cg_ids = sorted([
-                        bead.num + 1 for bead in mapping[mol]
-                        if bead.name in vbead.atoms
-                    ])
+                    cg_ids = sorted(
+                        [
+                            bead.num + 1
+                            for bead in mapping[mol]
+                            if bead.name in vbead.atoms
+                        ]
+                    )
                     cg_ids_string = " ".join(map(str, cg_ids))
-                    ret_lines.append("{0:^6d} {1:^6d} {2}".format(
-                        vbead.num + 1, vbead.gromacs_type_id, cg_ids_string))
-                    vsite_exclusions = "{} ".format(vbead.num +
-                                                    1) + cg_ids_string
+                    ret_lines.append(
+                        "{0:^6d} {1:^6d} {2}".format(
+                            vbead.num + 1, vbead.gromacs_type_id, cg_ids_string
+                        )
+                    )
+                    vsite_exclusions = "{} ".format(vbead.num + 1) + cg_ids_string
                     excl_lines.append(vsite_exclusions)
 
                 ret_lines.extend(excl_lines)
 
+            ret_lines.extend(write_bond_angle_dih(self.get_bond_lengths(mol), "bonds"))
             ret_lines.extend(
-                write_bond_angle_dih(self.get_bond_lengths(mol), "bonds"))
+                write_bond_angle_dih(self.get_bond_angles(mol), "angles", rad2deg=True)
+            )
             ret_lines.extend(
-                write_bond_angle_dih(self.get_bond_angles(mol),
-                                     "angles",
-                                     rad2deg=True))
+                write_bond_angle_dih(
+                    self.get_bond_dihedrals(mol),
+                    "dihedrals",
+                    multiplicity=1,
+                    rad2deg=True,
+                )
+            )
             ret_lines.extend(
-                write_bond_angle_dih(self.get_bond_dihedrals(mol),
-                                     "dihedrals",
-                                     multiplicity=1,
-                                     rad2deg=True))
-            ret_lines.extend(
-                write_bond_angle_dih(self.get_bond_length_constraints(mol),
-                                     "constraints",
-                                     print_fconst=False))
+                write_bond_angle_dih(
+                    self.get_bond_length_constraints(mol),
+                    "constraints",
+                    print_fconst=False,
+                )
+            )
 
         return ret_lines
 
@@ -461,25 +501,27 @@ class BondSet:
         calc = {
             2: mdtraj.compute_distances,
             3: mdtraj.compute_angles,
-            4: mdtraj.compute_dihedrals
+            4: mdtraj.compute_dihedrals,
         }
 
         for mol_name, mol_bonds in self._molecules.items():
             for bond in mol_bonds:
                 bond_indices = []
 
-                for prev_residue, residue, next_residue in sliding(
-                        frame.residues):
+                for prev_residue, residue, next_residue in sliding(frame.residues):
                     adj_res = {"-": prev_residue, "+": next_residue}
 
                     # Need access to adjacent residues, so can't just iterate over these directly
                     if residue.name == mol_name:
                         try:
-                            bond_indices.append([
-                                adj_res.get(atom_name[0], residue).atom(
-                                    atom_name.lstrip('-+')).index
-                                for atom_name in bond.atoms
-                            ])
+                            bond_indices.append(
+                                [
+                                    adj_res.get(atom_name[0], residue)
+                                    .atom(atom_name.lstrip("-+"))
+                                    .index
+                                    for atom_name in bond.atoms
+                                ]
+                            )
 
                         except AttributeError:
                             # AttributeError is raised when residues on end of chain calc bond to next
@@ -497,9 +539,9 @@ class BondSet:
                 pass
 
     @staticmethod
-    def _get_lines_for_bond_dump(bonds,
-                                 target_number: typing.Optional[int] = None,
-                                 rad2deg: bool = False) -> typing.List[str]:
+    def _get_lines_for_bond_dump(
+        bonds, target_number: typing.Optional[int] = None, rad2deg: bool = False
+    ) -> typing.List[str]:
         """Return a dump of bond measurements as a list of lines of text.
 
         :param bonds: Iterable of bonds
@@ -508,49 +550,53 @@ class BondSet:
         :return List[str]: Lines of bond measurements dump
         """
         ret_lines = []
-        for row in transpose_and_sample((bond.values for bond in bonds),
-                                        n=target_number):
+        for row in transpose_and_sample(
+            (bond.values for bond in bonds), n=target_number
+        ):
             if rad2deg:
                 row = [math.degrees(val) for val in row]
 
-            ret_lines.append((len(row) * '{:12.5f}').format(*row))
+            ret_lines.append((len(row) * "{:12.5f}").format(*row))
         return ret_lines
 
-    def dump_values(self, target_number: int = 10000,
-                    output_dir: typing.Optional[typing.Union[pathlib.Path, str]] = None) -> None:
+    def dump_values(
+        self,
+        target_number: int = 10000,
+        output_dir: typing.Optional[typing.Union[pathlib.Path, str]] = None,
+    ) -> None:
         """Output measured bond values to files for length, angles and dihedrals.
 
         :param int target_number: Approx number of sample measurements to output. If None, all samples will be output.
         :param output_dir: Directory to write output files to.
         """
         if output_dir is None:
-            output_dir = '.'
+            output_dir = "."
 
         # Cast to Path type
         output_dir = pathlib.Path(output_dir)
 
         for mol in self._molecules:
-            if mol == 'SOL':
+            if mol == "SOL":
                 continue
 
             bonds = self.get_bond_lengths(mol, with_constr=True)
             if bonds:
                 lines = BondSet._get_lines_for_bond_dump(bonds, target_number)
-                file_write_lines(output_dir.joinpath(f'{mol}_length.dat'), lines)
+                file_write_lines(output_dir.joinpath(f"{mol}_length.dat"), lines)
 
             bonds = self.get_bond_angles(mol)
             if bonds:
-                lines = BondSet._get_lines_for_bond_dump(bonds,
-                                                         target_number,
-                                                         rad2deg=True)
-                file_write_lines(output_dir.joinpath(f'{mol}_angle.dat'), lines)
+                lines = BondSet._get_lines_for_bond_dump(
+                    bonds, target_number, rad2deg=True
+                )
+                file_write_lines(output_dir.joinpath(f"{mol}_angle.dat"), lines)
 
             bonds = self.get_bond_dihedrals(mol)
             if bonds:
-                lines = BondSet._get_lines_for_bond_dump(bonds,
-                                                         target_number,
-                                                         rad2deg=True)
-                file_write_lines(output_dir.joinpath(f'{mol}_dihedral.dat'), lines)
+                lines = BondSet._get_lines_for_bond_dump(
+                    bonds, target_number, rad2deg=True
+                )
+                file_write_lines(output_dir.joinpath(f"{mol}_dihedral.dat"), lines)
 
     def __len__(self):
         return len(self._molecules)
diff --git a/pycgtool/forcefield.py b/pycgtool/forcefield.py
index be34b5727fe000f62e2a0b13ec82c766bfe0c5d6..8fb63398c1b2e6c7ae00791d675b27684a3fef0e 100644
--- a/pycgtool/forcefield.py
+++ b/pycgtool/forcefield.py
@@ -13,7 +13,9 @@ from .util import any_starts_with, backup_file, file_write_lines
 PathLike = typing.Union[pathlib.Path, str]
 
 
-def copy_files(src_dir: pathlib.Path, dest_dir: pathlib.Path, files: typing.Iterable[str]) -> None:
+def copy_files(
+    src_dir: pathlib.Path, dest_dir: pathlib.Path, files: typing.Iterable[str]
+) -> None:
     """Copy files from one directory to another."""
     for f in files:
         src_path = src_dir.joinpath(f)
@@ -25,37 +27,44 @@ class ForceField:
     """
     Class used to output a GROMACS .ff forcefield
     """
-    def __init__(self, name: str, dir_path: PathLike = pathlib.Path('.')):
+
+    def __init__(self, name: str, dir_path: PathLike = pathlib.Path(".")):
         """
         Open a named forcefield directory.  If it does not exist it is created.
 
         :param str name: Forcefield name to open/create
         """
-        self.directory = pathlib.Path(dir_path).joinpath(f'ff{name}.ff')
+        self.directory = pathlib.Path(dir_path).joinpath(f"ff{name}.ff")
         backup_file(self.directory)
         self.directory.mkdir(parents=True, exist_ok=True)
 
-        with open(self.directory.joinpath('forcefield.itp'), 'w') as itp:
-            print(f'#define _FF_PYCGTOOL_{name}', file=itp)
+        with open(self.directory.joinpath("forcefield.itp"), "w") as itp:
+            print(f"#define _FF_PYCGTOOL_{name}", file=itp)
             print('#include "martini_v2.2.itp"', file=itp)
 
-        data_dir = pathlib.Path(__file__).parent.joinpath('data')
+        data_dir = pathlib.Path(__file__).parent.joinpath("data")
 
         # Copy MARTINI files
-        copy_files(data_dir, self.directory, [
-            'martini_v2.2.itp',
-            'watermodels.dat',
-            'w.itp',
-        ])
+        copy_files(
+            data_dir,
+            self.directory,
+            [
+                "martini_v2.2.itp",
+                "watermodels.dat",
+                "w.itp",
+            ],
+        )
 
         # Create atomtypes.atp required for correct masses with pdb2gmx
-        atomtypes_atp = os.path.join(self.directory, 'atomtypes.atp')
-        with CFG(data_dir.joinpath('martini_v2.2.itp')) as itp, open(atomtypes_atp, 'w') as atomtypes:
-            for toks in itp['atomtypes']:
-                print(' '.join(toks), file=atomtypes)
+        atomtypes_atp = os.path.join(self.directory, "atomtypes.atp")
+        with CFG(data_dir.joinpath("martini_v2.2.itp")) as itp, open(
+            atomtypes_atp, "w"
+        ) as atomtypes:
+            for toks in itp["atomtypes"]:
+                print(" ".join(toks), file=atomtypes)
 
-        with open(self.directory.joinpath('forcefield.doc'), 'w') as doc:
-            print(f'PyCGTOOL produced MARTINI force field - {name}', file=doc)
+        with open(self.directory.joinpath("forcefield.doc"), "w") as doc:
+            print(f"PyCGTOOL produced MARTINI force field - {name}", file=doc)
 
     def write(self, filename: str, mapping, bonds):
         """
@@ -66,10 +75,10 @@ class ForceField:
         :param Iterable[Bond] bonds: CG Bonds object
         """
         lines, nterms, cterms = ForceField.write_rtp(mapping, bonds)
-        file_write_lines(self.directory.joinpath(f'{filename}.rtp'), lines)
+        file_write_lines(self.directory.joinpath(f"{filename}.rtp"), lines)
 
         lines = ForceField.write_r2b(nterms, cterms)
-        file_write_lines(self.directory.joinpath(f'{filename}.r2b'), lines)
+        file_write_lines(self.directory.joinpath(f"{filename}.r2b"), lines)
 
     @staticmethod
     def bond_section(bonds, section_header, multiplicity=None):
@@ -113,36 +122,38 @@ class ForceField:
         """
 
         def write_residue(mol_name, mol_mapping, strip=None, prepend=""):
-            ret_lines = [
-                "[ {0} ]".format(prepend + mol_name),
-                "  [ atoms ]"
-            ]
+            ret_lines = ["[ {0} ]".format(prepend + mol_name), "  [ atoms ]"]
 
             for bead in mol_mapping:
                 #                      name   type  charge  chg-group
-                ret_lines.append("    {:>4s} {:>4s} {:3.6f} {:4d}".format(
-                    bead.name, bead.type, bead.charge, 0
-                ))
-
-            for natoms, (section, multiplicity) in enumerate((("bonds", None),
-                                                              ("angles", None),
-                                                              ("dihedrals", 1)), start=2):
+                ret_lines.append(
+                    "    {:>4s} {:>4s} {:3.6f} {:4d}".format(
+                        bead.name, bead.type, bead.charge, 0
+                    )
+                )
+
+            for natoms, (section, multiplicity) in enumerate(
+                (("bonds", None), ("angles", None), ("dihedrals", 1)), start=2
+            ):
                 if strip is None:
                     bond_list = bonds.get_bonds(mol_name, natoms)
                 else:
-                    bond_list = bonds.get_bonds(mol_name, natoms, select=lambda bond: not any_starts_with(bond, strip))
+                    bond_list = bonds.get_bonds(
+                        mol_name,
+                        natoms,
+                        select=lambda bond: not any_starts_with(bond, strip),
+                    )
 
-                ret_lines.extend(ForceField.bond_section(bond_list, section, multiplicity))
+                ret_lines.extend(
+                    ForceField.bond_section(bond_list, section, multiplicity)
+                )
 
             return ret_lines
 
         n_terms = set()
         c_terms = set()
 
-        rtp_lines = [
-            "[ bondedtypes ]",
-            ("{:4d}" * 8).format(1, 1, 1, 1, 1, 1, 0, 0)
-        ]
+        rtp_lines = ["[ bondedtypes ]", ("{:4d}" * 8).format(1, 1, 1, 1, 1, 1, 0, 0)]
 
         for mol_name, mol_mapping in mapping.items():
             try:
@@ -153,14 +164,22 @@ class ForceField:
             needs_terminal_entry = ForceField.needs_terminal_entries(mol_name, bonds)
 
             if needs_terminal_entry[0]:
-                rtp_lines.extend(write_residue(mol_name, mol_mapping, strip="-", prepend="N"))
+                rtp_lines.extend(
+                    write_residue(mol_name, mol_mapping, strip="-", prepend="N")
+                )
                 n_terms.add(mol_name)
 
             if needs_terminal_entry[1]:
-                rtp_lines.extend(write_residue(mol_name, mol_mapping, strip="+", prepend="C"))
+                rtp_lines.extend(
+                    write_residue(mol_name, mol_mapping, strip="+", prepend="C")
+                )
                 c_terms.add(mol_name)
                 if needs_terminal_entry[0]:
-                    rtp_lines.extend(write_residue(mol_name, mol_mapping, strip=("-", "+"), prepend="2"))
+                    rtp_lines.extend(
+                        write_residue(
+                            mol_name, mol_mapping, strip=("-", "+"), prepend="2"
+                        )
+                    )
 
         return rtp_lines, n_terms, c_terms
 
@@ -177,13 +196,17 @@ class ForceField:
         """
         ret_lines = [
             "; rtp residue to rtp building block table",
-            ";     main  N-ter C-ter 2-ter"
+            ";     main  N-ter C-ter 2-ter",
         ]
 
         for resname in sorted(n_terms | c_terms):
             nter_str = ("N" + resname) if resname in n_terms else "-"
             cter_str = ("C" + resname) if resname in c_terms else "-"
             both_ter_str = ("2" + resname) if resname in (n_terms & c_terms) else "-"
-            ret_lines.append("{0:5s} {0:5s} {1:5s} {2:5s} {3:5s}".format(resname, nter_str, cter_str, both_ter_str))
+            ret_lines.append(
+                "{0:5s} {0:5s} {1:5s} {2:5s} {3:5s}".format(
+                    resname, nter_str, cter_str, both_ter_str
+                )
+            )
 
         return ret_lines
diff --git a/pycgtool/frame.py b/pycgtool/frame.py
index f1f303a706663cd50d06dfc3e71ce89f86674976..b4f7ae2320d5496d9059a897bab935b8b5cc393f 100644
--- a/pycgtool/frame.py
+++ b/pycgtool/frame.py
@@ -22,6 +22,7 @@ PathLike = typing.Union[pathlib.Path, str]
 
 class UnsupportedFormatException(Exception):
     """Exception raised when a topology/trajectory format cannot be parsed."""
+
     def __init__(self, msg=None):
         if msg is None:
             msg = "Topology/Trajectory format not supported by this reader"
@@ -30,6 +31,7 @@ class UnsupportedFormatException(Exception):
 
 class NonMatchingSystemError(ValueError):
     """Exception raised when topology and trajectory do not match."""
+
     def __init__(self, msg=None):
         if msg is None:
             msg = "Number of atoms does not match between topology and trajectory files"
@@ -40,9 +42,9 @@ def load_traj(filepath: PathLike, **kwargs) -> mdtraj.Trajectory:
     """Load a trajectory, if a PDB fails with zero box volume disable volume check and try again."""
     filepath = pathlib.Path(filepath)
 
-    if filepath.suffix.lower() == '.pdb':
+    if filepath.suffix.lower() == ".pdb":
         # PDBs have a couple of things that can go wrong - we handle these here...
-        kwargs.pop('top', None)  # Can't specify a topology for `load_pdb`
+        kwargs.pop("top", None)  # Can't specify a topology for `load_pdb`
 
         try:
             return mdtraj.load_pdb(str(filepath), **kwargs)
@@ -52,8 +54,8 @@ def load_traj(filepath: PathLike, **kwargs) -> mdtraj.Trajectory:
             # This can fail if the box volume is zero
             trajectory = mdtraj.load_pdb(str(filepath), no_boxchk=True, **kwargs)
             logger.warning(
-                'Unitcell has zero volume - periodic boundaries will not be accounted for. '
-                'If the molecule is split by a periodic boundary, results will be incorrect.'
+                "Unitcell has zero volume - periodic boundaries will not be accounted for. "
+                "If the molecule is split by a periodic boundary, results will be incorrect."
             )
 
             return trajectory
@@ -63,11 +65,14 @@ def load_traj(filepath: PathLike, **kwargs) -> mdtraj.Trajectory:
 
 class Frame:
     """Load and store data from a simulation trajectory."""
-    def __init__(self,
-                 topology_file: typing.Optional[PathLike] = None,
-                 trajectory_file: typing.Optional[PathLike] = None,
-                 frame_start: int = 0,
-                 frame_end: typing.Optional[int] = None):
+
+    def __init__(
+        self,
+        topology_file: typing.Optional[PathLike] = None,
+        trajectory_file: typing.Optional[PathLike] = None,
+        frame_start: int = 0,
+        frame_end: typing.Optional[int] = None,
+    ):
         """Load a simulation trajectory.
 
         :param topology_file: File containing system topology
@@ -77,19 +82,24 @@ class Frame:
         """
         if topology_file is not None:
             try:
-                logging.info('Loading topology file')
+                logging.info("Loading topology file")
                 self._trajectory = load_traj(topology_file)
                 self._topology = self._trajectory.topology
-                logging.info('Finished loading topology file')
+                logging.info("Finished loading topology file")
 
                 if trajectory_file is not None:
                     try:
-                        logging.info('Loading trajectory file - this may take a while')
-                        self._trajectory = load_traj(trajectory_file, top=self._topology)
-                        self._trajectory = self._slice_trajectory(frame_start, frame_end)
+                        logging.info("Loading trajectory file - this may take a while")
+                        self._trajectory = load_traj(
+                            trajectory_file, top=self._topology
+                        )
+                        self._trajectory = self._slice_trajectory(
+                            frame_start, frame_end
+                        )
                         logging.info(
-                            'Finished loading trajectory file - loaded %d frames',
-                            self._trajectory.n_frames)
+                            "Finished loading trajectory file - loaded %d frames",
+                            self._trajectory.n_frames,
+                        )
 
                     except ValueError as exc:
                         raise NonMatchingSystemError from exc
@@ -97,7 +107,7 @@ class Frame:
                 self._load_trajectory()
 
             except OSError as exc:
-                if 'no loader' in str(exc) or 'format is not supported' in str(exc):
+                if "no loader" in str(exc) or "format is not supported" in str(exc):
                     raise UnsupportedFormatException from exc
 
                 raise
@@ -107,9 +117,8 @@ class Frame:
             self._topology = mdtraj.Topology()
 
     def _slice_trajectory(
-            self,
-            frame_start: int = 0,
-            frame_end: typing.Optional[int] = None) -> mdtraj.Trajectory:
+        self, frame_start: int = 0, frame_end: typing.Optional[int] = None
+    ) -> mdtraj.Trajectory:
         """Slice input simulation trajectory to a subset of frames.
 
         :param frame_start: First frame of trajectory to use
@@ -177,20 +186,20 @@ class Frame:
         """Number of frames in the trajectory."""
         return self._trajectory.n_frames
 
-    def add_residue(self,
-                    name,
-                    chain: typing.Optional[
-                        mdtraj.core.topology.Residue] = None,
-                    **kwargs) -> mdtraj.core.topology.Residue:
+    def add_residue(
+        self,
+        name,
+        chain: typing.Optional[mdtraj.core.topology.Residue] = None,
+        **kwargs
+    ) -> mdtraj.core.topology.Residue:
         """Add a residue to the frame topology.
 
         :param name: Name of residue
         :param chain: MDTraj chain of residue
         :return: New residue
         """
-        if hasattr(self, '_trajectory'):
-            raise TypeError(
-                'Cannot edit residues if a trajectory has been loaded')
+        if hasattr(self, "_trajectory"):
+            raise TypeError("Cannot edit residues if a trajectory has been loaded")
 
         if chain is None:
             try:
@@ -202,8 +211,10 @@ class Frame:
         return self._topology.add_residue(name, chain, **kwargs)
 
     def add_atom(
-            self, name: str, element: typing.Optional[mdtraj.element.Element],
-            residue: mdtraj.core.topology.Residue
+        self,
+        name: str,
+        element: typing.Optional[mdtraj.element.Element],
+        residue: mdtraj.core.topology.Residue,
     ) -> mdtraj.core.topology.Atom:
         """Add an atom or CG bead to the frame topology.
 
@@ -212,16 +223,14 @@ class Frame:
         :param residue: MDTraj residue of atom
         :return: New atom
         """
-        if hasattr(self, '_trajectory'):
-            raise TypeError(
-                'Cannot edit atoms if a trajectory has been loaded')
+        if hasattr(self, "_trajectory"):
+            raise TypeError("Cannot edit atoms if a trajectory has been loaded")
 
         return self._topology.add_atom(name, element, residue)
 
-    def save(self,
-             filename: PathLike,
-             frame_number: typing.Optional[int] = None,
-             **kwargs) -> None:
+    def save(
+        self, filename: PathLike, frame_number: typing.Optional[int] = None, **kwargs
+    ) -> None:
         """Write trajctory to file.
 
         :param filename: Name of output file
@@ -249,14 +258,14 @@ class Frame:
 
         optional_values = {
             attr: getattr(self, attr, None)
-            for attr in {'time', 'unitcell_lengths', 'unitcell_angles'}
-        }  # yapf: disable
+            for attr in {"time", "unitcell_lengths", "unitcell_angles"}
+        }
 
-        new_trajectory = mdtraj.Trajectory(xyz,
-                                           topology=self._topology,
-                                           **optional_values)
+        new_trajectory = mdtraj.Trajectory(
+            xyz, topology=self._topology, **optional_values
+        )
 
-        if hasattr(self, '_trajectory'):
+        if hasattr(self, "_trajectory"):
             self._trajectory += new_trajectory
 
         else:
diff --git a/pycgtool/functionalforms.py b/pycgtool/functionalforms.py
index e55df800f8e026310ee4b8ea164f55037e288791..f237fe94a9c500620ddeec40faee3986a6e7a9d9 100644
--- a/pycgtool/functionalforms.py
+++ b/pycgtool/functionalforms.py
@@ -21,14 +21,17 @@ def get_functional_forms() -> enum.Enum:
         name = subclass.__name__
         enum_dict[name] = subclass
 
-    return enum.Enum('FunctionalForms', enum_dict)
+    return enum.Enum("FunctionalForms", enum_dict)
 
 
 class FunctionalForm(metaclass=abc.ABCMeta):
     """Parent class of any functional form used in Boltzmann Inversion to convert variance to a force constant."""
-    def __init__(self,
-                 mean_func: typing.Callable = np.nanmean,
-                 variance_func: typing.Callable = np.nanvar):
+
+    def __init__(
+        self,
+        mean_func: typing.Callable = np.nanmean,
+        variance_func: typing.Callable = np.nanvar,
+    ):
         """Inject functions for calculating the mean and variance into the Boltzmann Inversion equations.
 
         :param mean_func: Function to calculate the mean - default is np.nanmean
@@ -37,7 +40,9 @@ class FunctionalForm(metaclass=abc.ABCMeta):
         self._mean_func = mean_func
         self._variance_func = variance_func
 
-    def eqm(self, values: typing.Iterable[float], temp: float) -> float:  # pylint: disable=unused-argument
+    def eqm(
+        self, values: typing.Iterable[float], temp: float
+    ) -> float:  # pylint: disable=unused-argument
         """Calculate equilibrium value.
 
         May be overridden by functional forms.
@@ -78,7 +83,8 @@ class FunctionalForm(metaclass=abc.ABCMeta):
         if tipe is None:
             raise TypeError(
                 f"The functional form {cls.__name__} does not have a defined GROMACS "
-                f"potential type when used with {natoms} atoms.")
+                f"potential type when used with {natoms} atoms."
+            )
 
         return tipe
 
@@ -88,11 +94,12 @@ class Harmonic(FunctionalForm):
 
     See http://manual.gromacs.org/documentation/current/reference-manual/functions/bonded-interactions.html#harmonic-potential  # noqa
     """
+
     # TODO: Consider whether to use improper (type 2) instead, it is actually harmonic
     gromacs_type_ids = (1, 1, 1)
 
     def fconst(self, values: typing.Iterable[float], temp: float) -> float:
-        rt = 8.314 * temp / 1000.  # pylint: disable=invalid-name
+        rt = 8.314 * temp / 1000.0  # pylint: disable=invalid-name
         var = self._variance_func(values)
         return rt / var
 
@@ -104,34 +111,38 @@ class CosHarmonic(FunctionalForm):
 
     Uses the transformation in eqn 20 of the above source.
     """
+
     gromacs_type_ids = (None, 2, None)
 
     def fconst(self, values: typing.Iterable[float], temp: float) -> float:
-        rt = 8.314 * temp / 1000.  # pylint: disable=invalid-name
+        rt = 8.314 * temp / 1000.0  # pylint: disable=invalid-name
         mean = self.eqm(values, temp)
         var = self._variance_func(values)
-        return rt / (math.sin(mean)**2 * var)
+        return rt / (math.sin(mean) ** 2 * var)
 
 
 class MartiniDefaultLength(FunctionalForm):
     """Dummy functional form which returns a fixed force constant."""
+
     gromacs_type_ids = (1, None, None)
 
     def fconst(self, values: typing.Iterable[float], temp: float) -> float:
-        return 1250.
+        return 1250.0
 
 
 class MartiniDefaultAngle(FunctionalForm):
     """Dummy functional form which returns a fixed force constant."""
+
     gromacs_type_ids = (None, 2, None)
 
     def fconst(self, values: typing.Iterable[float], temp: float) -> float:
-        return 25.
+        return 25.0
 
 
 class MartiniDefaultDihedral(FunctionalForm):
     """Dummy functional form which returns a fixed force constant."""
+
     gromacs_type_ids = (None, None, 1)
 
     def fconst(self, values: typing.Iterable[float], temp: float) -> float:
-        return 50.
+        return 50.0
diff --git a/pycgtool/mapping.py b/pycgtool/mapping.py
index 03fd2ab3e68750c48fd6eb2f497210af5a3c780f..a02d3e3ef6b5248c20ff3ba5dfe7d0ae0fb27430 100644
--- a/pycgtool/mapping.py
+++ b/pycgtool/mapping.py
@@ -25,13 +25,16 @@ PathLike = typing.Union[str, pathlib.Path]
 
 class BeadMap:
     """Class holding values relating to the AA->CG transformation for a single bead."""
-    def __init__(self,
-                 name: str,
-                 num: int,
-                 type: str = None,
-                 atoms=None,
-                 charge: float = 0,
-                 mass: float = 0):
+
+    def __init__(
+        self,
+        name: str,
+        num: int,
+        type: str = None,
+        atoms=None,
+        charge: float = 0,
+        mass: float = 0,
+    ):
         """Create a single bead mapping.
 
         :param str name: The name of the bead
@@ -50,10 +53,8 @@ class BeadMap:
         self.atoms = atoms
         # NB: Mass weights are added in Mapping.__init__ if an itp file is provided
         self.weights_dict = {
-            "geom": np.array([1. / len(atoms) for _ in atoms],
-                             dtype=np.float32),
-            "first": np.array([1.] + [0. for _ in atoms[1:]],
-                              dtype=np.float32)
+            "geom": np.array([1.0 / len(atoms) for _ in atoms], dtype=np.float32),
+            "first": np.array([1.0] + [0.0 for _ in atoms[1:]], dtype=np.float32),
         }
         self.weights = self.weights_dict["geom"]
 
@@ -88,7 +89,8 @@ class BeadMap:
                 except KeyError:
                     raise RuntimeError(
                         f"Mass of atom {atom} could not be automatically assigned, "
-                        "map_center=mass is not available.")
+                        "map_center=mass is not available."
+                    )
 
             mass_array[i] = mass
 
@@ -97,7 +99,8 @@ class BeadMap:
         if not np.all(mass_array):
             raise RuntimeError(
                 "Some atom masses could not be automatically assigned, "
-                "map_center=mass is not available")
+                "map_center=mass is not available"
+            )
 
         mass_array /= self.mass
         self.weights_dict["mass"] = mass_array
@@ -113,12 +116,7 @@ class VirtualMap(BeadMap):
         :param List[str] atoms: The CG bead names from which the bead position is determined
         :param float charge: The net charge on the bead
         """
-        super().__init__(name,
-                         num,
-                         type=type,
-                         atoms=atoms,
-                         charge=charge,
-                         mass=0.)
+        super().__init__(name, num, type=type, atoms=atoms, charge=charge, mass=0.0)
 
         self.gromacs_type_id_dict = {"geom": 1, "mass": 2}
         self.gromacs_type_id = self.gromacs_type_id_dict["geom"]
@@ -133,7 +131,13 @@ class Mapping:
 
     Contains a dictionary of lists of BeadMaps.  Each list corresponds to a single molecule.
     """
-    def __init__(self, filename: PathLike, options, itp_filename: typing.Optional[PathLike] = None):
+
+    def __init__(
+        self,
+        filename: PathLike,
+        options,
+        itp_filename: typing.Optional[PathLike] = None,
+    ):
         """Read in the AA->CG mapping from a file.
 
         :param filename: File from which to read mapping
@@ -148,8 +152,7 @@ class Mapping:
 
         with CFG(filename) as cfg:
             for mol_name, mol_section in cfg.items():
-                mol_map, manual_charges = self._mol_map_from_section(
-                    mol_section)
+                mol_map, manual_charges = self._mol_map_from_section(mol_section)
 
                 self._mappings[mol_name] = mol_map
                 self._manual_charges[mol_name] = manual_charges
@@ -157,8 +160,9 @@ class Mapping:
         if itp_filename is not None:
             self._load_itp(itp_filename)
 
-        if ((self._map_center == "mass" or self._virtual_map_center == "mass")
-                and not self._masses_are_set):
+        if (
+            self._map_center == "mass" or self._virtual_map_center == "mass"
+        ) and not self._masses_are_set:
             self._guess_atom_masses()
 
         self._set_bead_weights()
@@ -176,16 +180,20 @@ class Mapping:
                     manual_charges = self._manual_charges[mol_name]
                     if manual_charges:
                         logger.warning(
-                            'Charges assigned in mapping for molecule %s, ignoring itp charges',
-                            mol_name)
+                            "Charges assigned in mapping for molecule %s, ignoring itp charges",
+                            mol_name,
+                        )
 
-                    self._itp_section_into_mol_map(mol_map, itp_mol_entry,
-                                                   manual_charges)
+                    self._itp_section_into_mol_map(
+                        mol_map, itp_mol_entry, manual_charges
+                    )
 
                 except KeyError:
                     logger.warning(
                         "No itp information for molecule %s found in %s",
-                        mol_name, itp.filename)
+                        mol_name,
+                        itp.filename,
+                    )
 
             self._masses_are_set = True
 
@@ -196,14 +204,16 @@ class Mapping:
                 if isinstance(bmap, VirtualMap):
                     bmap.weights = bmap.weights_dict[self._virtual_map_center]
                     bmap.gromacs_type_id = bmap.gromacs_type_id_dict[
-                        self._virtual_map_center]
+                        self._virtual_map_center
+                    ]
 
                 else:
                     bmap.weights = bmap.weights_dict[self._map_center]
 
     @staticmethod
-    def _itp_section_into_mol_map(mol_map: typing.List[BeadMap], itp_mol_entry,
-                                  manual_charges: bool) -> None:
+    def _itp_section_into_mol_map(
+        mol_map: typing.List[BeadMap], itp_mol_entry, manual_charges: bool
+    ) -> None:
         atoms = {}
         for toks in itp_mol_entry["atoms"]:
             # Store charge and mass
@@ -222,38 +232,41 @@ class Mapping:
 
         for bead in mol_map:
             if isinstance(bead, VirtualMap):
-                mass_array = np.array([
-                    real_bead.mass
-                    for real_bead in mol_map if real_bead.name in bead
-                ])
+                mass_array = np.array(
+                    [real_bead.mass for real_bead in mol_map if real_bead.name in bead]
+                )
                 weights_array = mass_array / sum(mass_array)
                 bead.weights_dict["mass"] = weights_array
 
                 if not manual_charges:
-                    bead.charge = sum([
-                        real_bead.charge for real_bead in mol_map
-                        if real_bead.name in bead
-                    ])
+                    bead.charge = sum(
+                        [
+                            real_bead.charge
+                            for real_bead in mol_map
+                            if real_bead.name in bead
+                        ]
+                    )
 
     @staticmethod
-    def _mol_map_from_section(
-            mol_section) -> typing.Tuple[typing.List[BeadMap], bool]:
+    def _mol_map_from_section(mol_section) -> typing.Tuple[typing.List[BeadMap], bool]:
         mol_map = []
         manual_charges = False
 
-        prefix_to_class = {'@v': VirtualMap}
+        prefix_to_class = {"@v": VirtualMap}
 
         for i, (name, typ, first, *atoms) in enumerate(mol_section):
             bead_class = BeadMap
             charge = 0
 
-            if name.startswith('@'):
+            if name.startswith("@"):
                 try:
                     bead_class = prefix_to_class[name]
                     name, typ, first, *atoms = mol_section[i][1:]
 
                 except KeyError as exc:
-                    raise ValueError(f'Invalid line prefix "{name}" in mapping file') from exc
+                    raise ValueError(
+                        f'Invalid line prefix "{name}" in mapping file'
+                    ) from exc
 
             try:
                 # Allow optional charge in mapping file
@@ -266,10 +279,9 @@ class Mapping:
                 atoms.insert(0, first)
 
             if not atoms:
-                raise ValueError(f'Bead {name} specification contains no atoms')
+                raise ValueError(f"Bead {name} specification contains no atoms")
 
-            mol_map.append(
-                bead_class(name, i, type=typ, atoms=atoms, charge=charge))
+            mol_map.append(bead_class(name, i, type=typ, atoms=atoms, charge=charge))
 
         return mol_map, manual_charges
 
@@ -300,8 +312,7 @@ class Mapping:
 
             for bead in mol_map:
                 bead.atoms = [
-                    atom_replacements[atom]
-                    if atom in atom_replacements else atom
+                    atom_replacements[atom] if atom in atom_replacements else atom
                     for atom in bead.atoms
                 ]
 
@@ -331,10 +342,14 @@ class Mapping:
             # Do this afterwards as it depends on real atom masses
             for bead in mol_mapping:
                 if isinstance(bead, VirtualMap):
-                    mass_array = np.array([
-                        real_bead.mass
-                        for real_bead in mol_mapping if real_bead.name in bead
-                    ], dtype=np.float32)  # yapf: disable
+                    mass_array = np.array(
+                        [
+                            real_bead.mass
+                            for real_bead in mol_mapping
+                            if real_bead.name in bead
+                        ],
+                        dtype=np.float32,
+                    )
 
                     weights_array = mass_array / sum(mass_array)
                     bead.weights_dict["mass"] = weights_array
@@ -342,12 +357,13 @@ class Mapping:
         self._masses_are_set = True
 
     def _cg_frame_setup(
-            self, aa_residues: typing.Iterable[mdtraj.core.topology.Residue]):
+        self, aa_residues: typing.Iterable[mdtraj.core.topology.Residue]
+    ):
         """Create a new CG Frame and populate beads
         :param aa_residues: Iterable of atomistic residues to map from
         :return: New CG Frame instance
         """
-        logger.info('Initialising output frame')
+        logger.info("Initialising output frame")
         cg_frame = Frame()
         missing_mappings = set()
 
@@ -360,7 +376,8 @@ class Mapping:
                     missing_mappings.add(aa_res.name)
                     logger.warning(
                         "A mapping has not been provided for '%s' residues, they will not be mapped",
-                        aa_res.name)
+                        aa_res.name,
+                    )
 
                 continue
 
@@ -368,7 +385,7 @@ class Mapping:
             for bmap in mol_map:
                 cg_frame.add_atom(bmap.name, None, cg_res)
 
-        logger.info('Finished initialising output frame')
+        logger.info("Finished initialising output frame")
         return cg_frame
 
     def apply(self, frame: Frame, cg_frame: typing.Optional[Frame] = None):
@@ -390,9 +407,8 @@ class Mapping:
         if not np.all(unitcell_lengths):
             unitcell_lengths = None
 
-        logger.info('Applying AA->CG mapping')
-        residues_to_map = (res for res in frame.residues
-                           if res.name in self._mappings)
+        logger.info("Applying AA->CG mapping")
+        residues_to_map = (res for res in frame.residues if res.name in self._mappings)
         for aa_res, cg_res in zip(residues_to_map, cg_frame.residues):
             mol_map = self._mappings[aa_res.name]
 
@@ -406,24 +422,24 @@ class Mapping:
 
                 else:
                     coords = np.array(
-                        [aa_res.atom(atom).coords for atom in bmap],
-                        dtype=np.float32)
-                    bead.coords = calc_coords_weight(ref_coords, coords,
-                                                     bmap.weights,
-                                                     unitcell_lengths)
+                        [aa_res.atom(atom).coords for atom in bmap], dtype=np.float32
+                    )
+                    bead.coords = calc_coords_weight(
+                        ref_coords, coords, bmap.weights, unitcell_lengths
+                    )
 
             for bead, bmap in zip(cg_res.atoms, mol_map):
                 if isinstance(bmap, VirtualMap):
                     coords = np.array(
-                        [cg_res.atom(atom).coords for atom in bmap],
-                        dtype=np.float32)
-                    bead.coords = calc_coords_weight(coords[0], coords,
-                                                     bmap.weights,
-                                                     unitcell_lengths)
+                        [cg_res.atom(atom).coords for atom in bmap], dtype=np.float32
+                    )
+                    bead.coords = calc_coords_weight(
+                        coords[0], coords, bmap.weights, unitcell_lengths
+                    )
 
         cg_frame.build_trajectory()
 
-        logger.info('Finished applying AA->CG mapping')
+        logger.info("Finished applying AA->CG mapping")
         return cg_frame
 
 
diff --git a/pycgtool/parsers/__init__.py b/pycgtool/parsers/__init__.py
index 99af7b3744a487be1861b5b4d8548c4be562fbcd..a454904887990bad117fe9d4db83377e665844d2 100644
--- a/pycgtool/parsers/__init__.py
+++ b/pycgtool/parsers/__init__.py
@@ -1,7 +1,4 @@
 from .cfg import CFG
 from .itp import ITP
 
-__all__ = [
-    'CFG',
-    'ITP'
-]
+__all__ = ["CFG", "ITP"]
diff --git a/pycgtool/parsers/cfg.py b/pycgtool/parsers/cfg.py
index 081fa697b415408003c40ee5063e35189553de06..060e1a63529acc161f6ba4043934d23b7d943e1b 100644
--- a/pycgtool/parsers/cfg.py
+++ b/pycgtool/parsers/cfg.py
@@ -13,6 +13,7 @@ PathLike = typing.Union[str, pathlib.Path]
 
 class DuplicateSectionError(KeyError):
     """Exception used to indicate that a section has appeared twice in a file."""
+
     def __init__(self, section, filename):
         msg = f"Section '{section}' appears twice in file '{filename}'."
         super().__init__(msg)
@@ -20,6 +21,7 @@ class DuplicateSectionError(KeyError):
 
 class NoSectionError(KeyError):
     """Exception used to indicate that a file contains no sections."""
+
     def __init__(self, filename):
         msg = f"File '{filename}' contains no '[]' section headers."
         super().__init__(msg)
@@ -30,6 +32,7 @@ class CFG(collections.OrderedDict, contextlib.AbstractContextManager):
 
     Contains a dictionary of Sections.
     """
+
     def __init__(self, filepath: typing.Optional[PathLike] = None):
         """Parse a config file and extract Sections."""
         super().__init__()
@@ -41,15 +44,15 @@ class CFG(collections.OrderedDict, contextlib.AbstractContextManager):
 
     def _read_line(self, line: str, filepath: pathlib.Path) -> str:
         # Strip comments
-        line = line.split(';')[0].strip()
+        line = line.split(";")[0].strip()
 
         # Handle include directive
-        if line.startswith('#include'):
+        if line.startswith("#include"):
             include_file = line.split()[1].strip('"')
             other = type(self)(filepath.parent.joinpath(include_file))
             self.update(other)
 
-            return ''  # Handle include then treat as empty line
+            return ""  # Handle include then treat as empty line
 
         return line
 
diff --git a/pycgtool/parsers/itp.py b/pycgtool/parsers/itp.py
index 6e19e7ddfea710783d6c4a3c165c1d1f1f053a2f..15b2a12bb6cc389b2c3fabf25317ad4e8d9440c3 100644
--- a/pycgtool/parsers/itp.py
+++ b/pycgtool/parsers/itp.py
@@ -14,6 +14,7 @@ class ITP(CFG):
 
     Contains a dictionary of sections for every molecule definition.
     """
+
     def _read_file(self, filepath: pathlib.Path) -> None:
         mol_sections = ["atoms", "bonds", "angles", "dihedrals"]
 
@@ -54,5 +55,8 @@ class ITP(CFG):
                     raise NoSectionError(filepath)
 
                 else:
-                    logger.info("File '%s' contains unexpected section '%s'",
-                                filepath, curr_section)
+                    logger.info(
+                        "File '%s' contains unexpected section '%s'",
+                        filepath,
+                        curr_section,
+                    )
diff --git a/pycgtool/util.py b/pycgtool/util.py
index c308ca95b0a984761964e9c35deeffc65e37db0c..7ec74bc3ddaf823c128920a42f867eac91a2d0d9 100644
--- a/pycgtool/util.py
+++ b/pycgtool/util.py
@@ -84,10 +84,16 @@ def extend_graph_chain(extend, pairs):
                 if node2.startswith("+"):
                     for pair2 in pairs:
                         if node2.strip("+") == pair2[0] and "+" + pair2[1] not in chain:
-                            append_if_not_in(ret, spare + (node1, node2, "+" + pair2[1]))
+                            append_if_not_in(
+                                ret, spare + (node1, node2, "+" + pair2[1])
+                            )
 
-                        elif node2.strip("+") == pair2[1] and "+" + pair2[0] not in chain:
-                            append_if_not_in(ret, spare + (node1, node2, "+" + pair2[0]))
+                        elif (
+                            node2.strip("+") == pair2[1] and "+" + pair2[0] not in chain
+                        ):
+                            append_if_not_in(
+                                ret, spare + (node1, node2, "+" + pair2[0])
+                            )
 
             except AttributeError:
                 pass
@@ -98,7 +104,9 @@ def extend_graph_chain(extend, pairs):
     return ret
 
 
-def transpose_and_sample(sequence: typing.Iterable, n: typing.Optional[int] = None) -> typing.List:
+def transpose_and_sample(
+    sequence: typing.Iterable, n: typing.Optional[int] = None
+) -> typing.List:
     """Transpose a sequence of lists and sample to provide target number of rows.
 
     Skip rows containing non-finite numbers (e.g. NaN).
@@ -106,6 +114,7 @@ def transpose_and_sample(sequence: typing.Iterable, n: typing.Optional[int] = No
     :param sequence: 2d sequence object to transpose
     :param n: Number of samples to take
     """
+
     def no_nans(row):
         return np.isfinite(row).all()
 
@@ -143,12 +152,12 @@ def backup_file(path: PathLike) -> pathlib.Path:
     counter = 1
 
     while new_path.exists():
-        new_path = path.parent.joinpath(f'#{path.name}.{counter}#')
+        new_path = path.parent.joinpath(f"#{path.name}.{counter}#")
         counter += 1
 
     if new_path != path:
         path.rename(new_path)
-        logger.warning('Existing file %s backed up as %s', path.name, new_path.name)
+        logger.warning("Existing file %s backed up as %s", path.name, new_path.name)
 
     return new_path
 
@@ -166,7 +175,7 @@ def sliding(vals: typing.Iterable):
         current = next(it)
 
     except StopIteration:
-        raise EmptyIterableError('Cannot make sliding window over empty iterable')
+        raise EmptyIterableError("Cannot make sliding window over empty iterable")
 
     for nxt in it:
         yield (prev, current, nxt)
@@ -221,7 +230,9 @@ def cmp_whitespace_float(ref_lines, test_lines, rtol=0.01, verbose=False):
     :return: True if all lines are the same, else False
     """
     diff_lines = []
-    for i, (ref_line, test_line) in enumerate(itertools.zip_longest(ref_lines, test_lines)):
+    for i, (ref_line, test_line) in enumerate(
+        itertools.zip_longest(ref_lines, test_lines)
+    ):
         # Shortcut trivial comparisons
         if ref_line is None or test_line is None:
             diff_lines.append((i, ref_line, test_line))
@@ -234,8 +245,9 @@ def cmp_whitespace_float(ref_lines, test_lines, rtol=0.01, verbose=False):
         test_toks = test_line.split()
 
         # Check for float comparison of tokens on line
-        for ref_tok, test_tok in itertools.zip_longest(map(number_or_string, ref_toks),
-                                                       map(number_or_string, test_toks)):
+        for ref_tok, test_tok in itertools.zip_longest(
+            map(number_or_string, ref_toks), map(number_or_string, test_toks)
+        ):
             if ref_tok != test_tok:
                 try:
                     if abs(ref_tok - test_tok) > abs(ref_tok) * rtol:
@@ -287,8 +299,9 @@ def any_starts_with(iterable: typing.Iterable, char: str) -> bool:
 
 
 def load_optional_topology(
-        trajfile: typing.Union[str, pathlib.Path],
-        topfile: typing.Optional[typing.Union[str, pathlib.Path]] = None) -> mdtraj.Trajectory:
+    trajfile: typing.Union[str, pathlib.Path],
+    topfile: typing.Optional[typing.Union[str, pathlib.Path]] = None,
+) -> mdtraj.Trajectory:
     """Load an MDTraj trajectory with or without a separate topology file.
 
     :param trajfile: Trajectory file
@@ -301,9 +314,11 @@ def load_optional_topology(
     return mdtraj.load(str(trajfile), top=str(topfile))
 
 
-def compare_trajectories(*trajectory_files: typing.Union[str, pathlib.Path],
-                         topology_file: typing.Optional[typing.Union[str, pathlib.Path]] = None,
-                         rtol: float = 0.001) -> bool:
+def compare_trajectories(
+    *trajectory_files: typing.Union[str, pathlib.Path],
+    topology_file: typing.Optional[typing.Union[str, pathlib.Path]] = None,
+    rtol: float = 0.001,
+) -> bool:
     """Compare multiple trajectory files for equality.
 
     :param trajectory_files: Paths of trajectory file to compare
@@ -317,22 +332,22 @@ def compare_trajectories(*trajectory_files: typing.Union[str, pathlib.Path],
             traj = load_optional_topology(traj_file, topology_file)
 
         except ValueError as exc:
-            raise ValueError('Topology does not match') from exc
+            raise ValueError("Topology does not match") from exc
 
         # Conditions from Trajectory.__hash__
         if traj.topology != ref_traj.topology:
-            raise ValueError('Topology does not match')
+            raise ValueError("Topology does not match")
 
         if len(traj) != len(ref_traj):
-            raise ValueError('Length does not match')
+            raise ValueError("Length does not match")
 
         if not np.allclose(traj.xyz, ref_traj.xyz, rtol=rtol):
-            raise ValueError('Coordinates do not match')
+            raise ValueError("Coordinates do not match")
 
         if not np.allclose(traj.time, ref_traj.time, rtol=rtol):
-            raise ValueError('Time does not match')
+            raise ValueError("Time does not match")
 
         if not np.allclose(traj.unitcell_vectors, ref_traj.unitcell_vectors, rtol=rtol):
-            raise ValueError('Unitcell does not match')
+            raise ValueError("Unitcell does not match")
 
     return True
diff --git a/pyproject.toml b/pyproject.toml
index ac464376ababd877acd90efb0ca2d279cc8a72b7..bd967815f08ac9f3d6f15f1714db6bf61d2622da 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -26,7 +26,7 @@ classifiers = [
 pycgtool = "pycgtool.__main__:main"
 
 [tool.poetry.dependencies]
-python = "^3.6"
+python = "^3.6.2"
 wheel = "^0.35.1"
 numpy = [
     { version = "^1.19.1", python = "<3.7" },  # Support for Python 3.6 was dropped in 1.20
@@ -48,7 +48,6 @@ myst-parser = { version = "^0.13.5", optional = true }
 prospector = "^1.3.0"
 pytest = "^6.0.1"
 pytest-cov = "^2.10.1"
-yapf = "^0.31.0"
 vulture = "^2.1"
 rstcheck = "^3.3.1"
 flake8 = "^3.8.4"
@@ -57,14 +56,12 @@ sphinx-autoapi = "^1.5.1"
 sphinx-rtd-theme = "^0.5.1"
 myst-parser = "^0.13.5"
 tox = "^3.23.0"
+black = "^22.1.0"
 
 [tool.poetry.extras]
 backmapping = ["mdplus"]
 docs = ["Sphinx", "sphinx-autoapi", "sphinx-rtd-theme", "myst-parser"]
 
-[tool.yapf]
-column_limit = 110
-
 [build-system]
 requires = ["poetry-core>=1.0.0"]
 build-backend = "poetry.core.masonry.api"
diff --git a/tests/test_bondset.py b/tests/test_bondset.py
index e98f742265903aff89563eeb8c5ac78012540c30..b4bf52cbbb5b8a610459ef5da7881bfb286a2a34 100644
--- a/tests/test_bondset.py
+++ b/tests/test_bondset.py
@@ -1,4 +1,3 @@
-
 import math
 import os
 import pathlib
@@ -13,8 +12,8 @@ from pycgtool.util import cmp_file_whitespace_float
 
 class DummyOptions:
     constr_threshold = 100000
-    map_center = 'geom'
-    virtual_map_center = 'geom'
+    map_center = "geom"
+    virtual_map_center = "geom"
     angle_default_fc = False
     generate_angles = True
     generate_dihedrals = False
@@ -34,9 +33,7 @@ class DummyBond:
 # TODO add setup and teardown to put all files in a tmp directory
 class BondSetTest(unittest.TestCase):
     base_dir = pathlib.Path(__file__).absolute().parent
-    data_dir = base_dir.joinpath('data')
-
-    # Columns are: eqm value, standard fc, MARTINI defaults fc
+    data_dir = base_dir.joinpath("data")
     invert_test_ref_data = [
         ( 0.220474419132,  72658.18163, 1250),
         ( 0.221844516963,  64300.01188, 1250),
@@ -55,40 +52,43 @@ class BondSetTest(unittest.TestCase):
         ( 2.768063749729,    135.50927,   50),
         (-2.213755550432,     51.13975,   50),
         ( 1.456495664734,     59.38162,   50),
-        (-1.826134298998,    279.80889,   50)
-    ]
+        (-1.826134298998,    279.80889,   50),
+    ]  # fmt: skip
 
     def test_bondset_create(self):
-        measure = BondSet(self.data_dir.joinpath('sugar.bnd'), DummyOptions)
+        measure = BondSet(self.data_dir.joinpath("sugar.bnd"), DummyOptions)
         self.assertEqual(1, len(measure))
-        self.assertTrue('ALLA' in measure)
-        self.assertEqual(18, len(measure['ALLA']))
+        self.assertTrue("ALLA" in measure)
+        self.assertEqual(18, len(measure["ALLA"]))
 
     def test_bondset_apply(self):
-        measure = BondSet(self.data_dir.joinpath('sugar.bnd'), DummyOptions)
-        frame = Frame(self.data_dir.joinpath('sugar-cg.gro'))
+        measure = BondSet(self.data_dir.joinpath("sugar.bnd"), DummyOptions)
+        frame = Frame(self.data_dir.joinpath("sugar-cg.gro"))
         measure.apply(frame)
 
         # First six are bond lengths
-        self.assertEqual(1, len(measure['ALLA'][0].values))
-        self.assertAlmostEqual(0.2225376, measure['ALLA'][0].values[0],
-                               delta=0.2225376 / 500)
+        self.assertEqual(1, len(measure["ALLA"][0].values))
+        self.assertAlmostEqual(
+            0.2225376, measure["ALLA"][0].values[0], delta=0.2225376 / 500
+        )
         # Second six are angles
-        self.assertEqual(1, len(measure['ALLA'][6].values))
+        self.assertEqual(1, len(measure["ALLA"][6].values))
         expected = math.radians(77.22779289)
-        self.assertAlmostEqual(expected, measure['ALLA'][6].values[0],
-                               delta=expected / 500)
+        self.assertAlmostEqual(
+            expected, measure["ALLA"][6].values[0], delta=expected / 500
+        )
         # Final six are dihedrals
-        self.assertEqual(1, len(measure['ALLA'][12].values))
+        self.assertEqual(1, len(measure["ALLA"][12].values))
         expected = math.radians(-89.5552903)
-        self.assertAlmostEqual(expected, measure['ALLA'][12].values[0],
-                               delta=abs(expected) / 500)
+        self.assertAlmostEqual(
+            expected, measure["ALLA"][12].values[0], delta=abs(expected) / 500
+        )
 
     def test_bondset_remove_triangles(self):
-        bondset = BondSet(self.data_dir.joinpath('triangle.bnd'), DummyOptions)
-        angles = bondset.get_bond_angles('TRI', exclude_triangle=False)
+        bondset = BondSet(self.data_dir.joinpath("triangle.bnd"), DummyOptions)
+        angles = bondset.get_bond_angles("TRI", exclude_triangle=False)
         self.assertEqual(3, len(angles))
-        angles = bondset.get_bond_angles('TRI', exclude_triangle=True)
+        angles = bondset.get_bond_angles("TRI", exclude_triangle=True)
         self.assertEqual(0, len(angles))
 
     def support_check_mean_fc(self, mol_bonds, fc_column_number):
@@ -98,87 +98,95 @@ class BondSetTest(unittest.TestCase):
 
         for i, bond in enumerate(mol_bonds):
             ref = self.invert_test_ref_data
-            self.assertAlmostEqual(ref[i][0], bond.eqm,
-                                   delta=abs(ref[i][0] * accuracy))
-            self.assertAlmostEqual(ref[i][fc_column_number], bond.fconst,
-                                   delta=abs(ref[i][fc_column_number] * accuracy))
+            self.assertAlmostEqual(ref[i][0], bond.eqm, delta=abs(ref[i][0] * accuracy))
+            self.assertAlmostEqual(
+                ref[i][fc_column_number],
+                bond.fconst,
+                delta=abs(ref[i][fc_column_number] * accuracy),
+            )
 
     def test_bondset_boltzmann_invert(self):
-        measure = BondSet(self.data_dir.joinpath('sugar.bnd'), DummyOptions)
-        frame = Frame(self.data_dir.joinpath('sugar.gro'), self.data_dir.joinpath('sugar.xtc'))
-        mapping = Mapping(self.data_dir.joinpath('sugar.map'), DummyOptions)
+        measure = BondSet(self.data_dir.joinpath("sugar.bnd"), DummyOptions)
+        frame = Frame(
+            self.data_dir.joinpath("sugar.gro"), self.data_dir.joinpath("sugar.xtc")
+        )
+        mapping = Mapping(self.data_dir.joinpath("sugar.map"), DummyOptions)
 
         cg_frame = mapping.apply(frame)
 
         measure.apply(cg_frame)
         measure.boltzmann_invert()
 
-        self.support_check_mean_fc(measure['ALLA'], 1)
+        self.support_check_mean_fc(measure["ALLA"], 1)
 
     def test_bondset_boltzmann_invert_default_fc(self):
         class DefaultOptions(DummyOptions):
             default_fc = True
 
-        measure = BondSet(self.data_dir.joinpath('sugar.bnd'), DefaultOptions)
-        frame = Frame(self.data_dir.joinpath('sugar.gro'), self.data_dir.joinpath('sugar.xtc'))
-        mapping = Mapping(self.data_dir.joinpath('sugar.map'), DefaultOptions)
+        measure = BondSet(self.data_dir.joinpath("sugar.bnd"), DefaultOptions)
+        frame = Frame(
+            self.data_dir.joinpath("sugar.gro"), self.data_dir.joinpath("sugar.xtc")
+        )
+        mapping = Mapping(self.data_dir.joinpath("sugar.map"), DefaultOptions)
 
         cg_frame = mapping.apply(frame)
 
         measure.apply(cg_frame)
         measure.boltzmann_invert()
 
-        self.support_check_mean_fc(measure['ALLA'], 2)
+        self.support_check_mean_fc(measure["ALLA"], 2)
 
     def test_bondset_boltzmann_invert_manual_default_fc(self):
         class FuncFormOptions(DummyOptions):
-            length_form = 'MartiniDefaultLength'
-            angle_form = 'MartiniDefaultAngle'
-            dihedral_form = 'MartiniDefaultDihedral'
+            length_form = "MartiniDefaultLength"
+            angle_form = "MartiniDefaultAngle"
+            dihedral_form = "MartiniDefaultDihedral"
 
-        measure = BondSet(self.data_dir.joinpath('sugar.bnd'), FuncFormOptions)
-        frame = Frame(self.data_dir.joinpath('sugar.gro'), self.data_dir.joinpath('sugar.xtc'))
-        mapping = Mapping(self.data_dir.joinpath('sugar.map'), DummyOptions)
+        measure = BondSet(self.data_dir.joinpath("sugar.bnd"), FuncFormOptions)
+        frame = Frame(
+            self.data_dir.joinpath("sugar.gro"), self.data_dir.joinpath("sugar.xtc")
+        )
+        mapping = Mapping(self.data_dir.joinpath("sugar.map"), DummyOptions)
 
         cg_frame = mapping.apply(frame)
 
         measure.apply(cg_frame)
         measure.boltzmann_invert()
 
-        self.support_check_mean_fc(measure['ALLA'], 2)
+        self.support_check_mean_fc(measure["ALLA"], 2)
 
     def test_bondset_polymer(self):
-        bondset = BondSet(self.data_dir.joinpath('polyethene.bnd'), DummyOptions)
-        frame = Frame(self.data_dir.joinpath('polyethene.gro'))
+        bondset = BondSet(self.data_dir.joinpath("polyethene.bnd"), DummyOptions)
+        frame = Frame(self.data_dir.joinpath("polyethene.gro"))
         bondset.apply(frame)
 
-        self.assertEqual(5, len(bondset['ETH'][0].values))
-        self.assertEqual(4, len(bondset['ETH'][1].values))
-        self.assertEqual(4, len(bondset['ETH'][2].values))
-        self.assertEqual(4, len(bondset['ETH'][3].values))
+        self.assertEqual(5, len(bondset["ETH"][0].values))
+        self.assertEqual(4, len(bondset["ETH"][1].values))
+        self.assertEqual(4, len(bondset["ETH"][2].values))
+        self.assertEqual(4, len(bondset["ETH"][3].values))
 
         bondset.boltzmann_invert()
 
-        self.assertAlmostEqual(0.107, bondset['ETH'][0].eqm,
-                               delta=0.107 / 500)
-        self.assertAlmostEqual(0.107, bondset['ETH'][1].eqm,
-                               delta=0.107 / 500)
+        self.assertAlmostEqual(0.107, bondset["ETH"][0].eqm, delta=0.107 / 500)
+        self.assertAlmostEqual(0.107, bondset["ETH"][1].eqm, delta=0.107 / 500)
 
     def test_bondset_pbc(self):
-        bondset = BondSet(self.data_dir.joinpath('polyethene.bnd'), DummyOptions)
-        frame = Frame(self.data_dir.joinpath('pbcpolyethene.gro'))
+        bondset = BondSet(self.data_dir.joinpath("polyethene.bnd"), DummyOptions)
+        frame = Frame(self.data_dir.joinpath("pbcpolyethene.gro"))
 
         bondset.apply(frame)
         bondset.boltzmann_invert()
 
-        for bond in bondset.get_bond_lengths('ETH', True):
-            self.assertAlmostEqual(1., bond.eqm)
+        for bond in bondset.get_bond_lengths("ETH", True):
+            self.assertAlmostEqual(1.0, bond.eqm)
             self.assertEqual(math.inf, bond.fconst)
 
     def test_full_itp_sugar(self):
-        measure = BondSet(self.data_dir.joinpath('sugar.bnd'), DummyOptions)
-        frame = Frame(self.data_dir.joinpath('sugar.gro'), self.data_dir.joinpath('sugar.xtc'))
-        mapping = Mapping(self.data_dir.joinpath('sugar.map'), DummyOptions)
+        measure = BondSet(self.data_dir.joinpath("sugar.bnd"), DummyOptions)
+        frame = Frame(
+            self.data_dir.joinpath("sugar.gro"), self.data_dir.joinpath("sugar.xtc")
+        )
+        mapping = Mapping(self.data_dir.joinpath("sugar.map"), DummyOptions)
 
         cg_frame = mapping.apply(frame)
 
@@ -188,23 +196,25 @@ class BondSetTest(unittest.TestCase):
         with tempfile.TemporaryDirectory() as tmpdir:
             tmp_path = pathlib.Path(tmpdir)
 
-            measure.write_itp(tmp_path.joinpath('sugar_out.itp'), mapping)
+            measure.write_itp(tmp_path.joinpath("sugar_out.itp"), mapping)
 
             self.assertTrue(
                 cmp_file_whitespace_float(
-                    tmp_path.joinpath('sugar_out.itp'),
-                    self.data_dir.joinpath('sugar_out.itp'),
+                    tmp_path.joinpath("sugar_out.itp"),
+                    self.data_dir.joinpath("sugar_out.itp"),
                     rtol=0.005,
-                    verbose=True))
+                    verbose=True,
+                )
+            )
 
     def test_full_itp_vsites(self):
         """Test full operation to output of .itp file for molecule with vsites."""
         options = DummyOptions()
         options.generate_angles = False
 
-        measure = BondSet(self.data_dir.joinpath('martini3/naphthalene.bnd'), options)
-        frame = Frame(self.data_dir.joinpath('martini3/naphthalene.gro'))
-        mapping = Mapping(self.data_dir.joinpath('martini3/naphthalene.map'), options)
+        measure = BondSet(self.data_dir.joinpath("martini3/naphthalene.bnd"), options)
+        frame = Frame(self.data_dir.joinpath("martini3/naphthalene.gro"))
+        mapping = Mapping(self.data_dir.joinpath("martini3/naphthalene.map"), options)
 
         cg_frame = mapping.apply(frame)
 
@@ -213,24 +223,27 @@ class BondSetTest(unittest.TestCase):
 
         with tempfile.TemporaryDirectory() as tmpdir:
             tmp_path = pathlib.Path(tmpdir)
-            measure.write_itp(tmp_path.joinpath('naphthalene_out.itp'),
-                              mapping)
+            measure.write_itp(tmp_path.joinpath("naphthalene_out.itp"), mapping)
 
             self.assertTrue(
                 cmp_file_whitespace_float(
-                    tmp_path.joinpath('naphthalene_out.itp'),
-                    self.data_dir.joinpath('martini3/naphthalene_out.itp'),
+                    tmp_path.joinpath("naphthalene_out.itp"),
+                    self.data_dir.joinpath("martini3/naphthalene_out.itp"),
                     rtol=0.005,
-                    verbose=True))
+                    verbose=True,
+                )
+            )
 
     def test_duplicate_atoms_in_bond(self):
         with self.assertRaises(ValueError):
-            _ = BondSet(self.data_dir.joinpath('duplicate_atoms.bnd'), DummyOptions)
+            _ = BondSet(self.data_dir.joinpath("duplicate_atoms.bnd"), DummyOptions)
 
     def test_dump_bonds(self):
-        measure = BondSet(self.data_dir.joinpath('sugar.bnd'), DummyOptions)
-        frame = Frame(self.data_dir.joinpath('sugar.gro'), self.data_dir.joinpath('sugar.xtc'))
-        mapping = Mapping(self.data_dir.joinpath('sugar.map'), DummyOptions)
+        measure = BondSet(self.data_dir.joinpath("sugar.bnd"), DummyOptions)
+        frame = Frame(
+            self.data_dir.joinpath("sugar.gro"), self.data_dir.joinpath("sugar.xtc")
+        )
+        mapping = Mapping(self.data_dir.joinpath("sugar.map"), DummyOptions)
 
         cg_frame = mapping.apply(frame)
 
@@ -238,27 +251,27 @@ class BondSetTest(unittest.TestCase):
         measure.boltzmann_invert()
         measure.dump_values()
 
-        filenames = ('ALLA_length.dat', 'ALLA_angle.dat', 'ALLA_dihedral.dat')
+        filenames = ("ALLA_length.dat", "ALLA_angle.dat", "ALLA_dihedral.dat")
         for filename in filenames:
             self.assertTrue(
-                cmp_file_whitespace_float(self.data_dir.joinpath(filename),
-                                          filename,
-                                          rtol=0.008,
-                                          verbose=True))
+                cmp_file_whitespace_float(
+                    self.data_dir.joinpath(filename), filename, rtol=0.008, verbose=True
+                )
+            )
             os.remove(filename)
 
     def test_get_lines_for_bond_dump(self):
         expected = [
-            '     0.00000     1.00000     2.00000',
-            '     1.00000     2.00000     3.00000',
-            '     2.00000     3.00000     4.00000',
-            '     3.00000     4.00000     5.00000'
+            "     0.00000     1.00000     2.00000",
+            "     1.00000     2.00000     3.00000",
+            "     2.00000     3.00000     4.00000",
+            "     3.00000     4.00000     5.00000",
         ]
 
         bonds = [
             DummyBond(None, None, None, values=[0, 1, 2, 3]),
             DummyBond(None, None, None, values=[1, 2, 3, 4]),
-            DummyBond(None, None, None, values=[2, 3, 4, 5])
+            DummyBond(None, None, None, values=[2, 3, 4, 5]),
         ]
 
         output = BondSet._get_lines_for_bond_dump(bonds)
@@ -267,16 +280,16 @@ class BondSetTest(unittest.TestCase):
 
     def test_get_lines_for_bond_dump_sample(self):
         expected = [
-            '     0.00000     1.00000     2.00000',
-            '     1.00000     2.00000     3.00000',
-            '     2.00000     3.00000     4.00000',
-            '     3.00000     4.00000     5.00000'
+            "     0.00000     1.00000     2.00000",
+            "     1.00000     2.00000     3.00000",
+            "     2.00000     3.00000     4.00000",
+            "     3.00000     4.00000     5.00000",
         ]
 
         bonds = [
             DummyBond(None, None, None, values=[0, 1, 2, 3]),
             DummyBond(None, None, None, values=[1, 2, 3, 4]),
-            DummyBond(None, None, None, values=[2, 3, 4, 5])
+            DummyBond(None, None, None, values=[2, 3, 4, 5]),
         ]
 
         nlines = 2
diff --git a/tests/test_forcefield.py b/tests/test_forcefield.py
index 57d5ed869f3d70e42f600c6627320f91158cbb56..30e0c11304976c9228972869bec0c8a6a72f6236 100644
--- a/tests/test_forcefield.py
+++ b/tests/test_forcefield.py
@@ -29,7 +29,9 @@ class DummyBondSet:
     def get_bonds(self, mol, natoms, select=lambda x: True):
         if natoms == -1:
             return [bond for bond in self.bonds if select(bond)]
-        return [bond for bond in self.bonds if len(bond.atoms) == natoms and select(bond)]
+        return [
+            bond for bond in self.bonds if len(bond.atoms) == natoms and select(bond)
+        ]
 
     def get_bond_lengths(self, *args, **kwargs):
         return self.get_bonds(None, 2)
@@ -44,16 +46,18 @@ class DummyBondSet:
 class ForceFieldTest(unittest.TestCase):
     def setUp(self):
         self.bonds = [
-            DummyBond(["a",  "b"], 1, 100),
-            DummyBond(["b",  "c"], 2,  50),
-            DummyBond(["c", "+a"], 3,  20),
+            DummyBond(["a", "b"], 1, 100),
+            DummyBond(["b", "c"], 2, 50),
+            DummyBond(["c", "+a"], 3, 20),
         ]
 
-        self.mapping = {"Dummy": [
-                DummyBMap("a", "a1",  1),
+        self.mapping = {
+            "Dummy": [
+                DummyBMap("a", "a1", 1),
                 DummyBMap("b", "b2", -1),
-                DummyBMap("c", "c3",  0),
-        ]}
+                DummyBMap("c", "c3", 0),
+            ]
+        }
 
         self.bondset = DummyBondSet(self.bonds, "Dummy")
 
@@ -62,7 +66,7 @@ class ForceFieldTest(unittest.TestCase):
             tmp_dir = pathlib.Path(t)
 
             name = "test"
-            ff_dir = tmp_dir.joinpath('fftest.ff')
+            ff_dir = tmp_dir.joinpath("fftest.ff")
 
             ForceField(name, dir_path=tmp_dir)
             self.assertTrue(ff_dir.exists())
@@ -76,7 +80,7 @@ class ForceFieldTest(unittest.TestCase):
             "  [ section ]",
             "       a    b      1.00000    100.00000",
             "       b    c      2.00000     50.00000",
-            "       c   +a      3.00000     20.00000"
+            "       c   +a      3.00000     20.00000",
         ]
 
         self.assertListEqual(expected, ForceField.bond_section(self.bonds, "section"))
@@ -91,7 +95,7 @@ class ForceFieldTest(unittest.TestCase):
             "a     a     Na    Ca    2a   ",
             "b     b     Nb    -     -    ",
             "c     c     -     Cc    -    ",
-            "d     d     Nd    Cd    2d   "
+            "d     d     Nd    Cd    2d   ",
         ]
 
         output = ForceField.write_r2b(nter, cter)
@@ -117,7 +121,7 @@ class ForceFieldTest(unittest.TestCase):
             "       c   c3 0.000000    0",
             "  [ bonds ]",
             "       a    b      1.00000    100.00000",
-            "       b    c      2.00000     50.00000"
+            "       b    c      2.00000     50.00000",
         ]
 
         output, nter, cter = ForceField.write_rtp(self.mapping, self.bondset)
@@ -135,5 +139,5 @@ class ForceFieldTest(unittest.TestCase):
         self.assertTrue(cter)
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     unittest.main()
diff --git a/tests/test_frame.py b/tests/test_frame.py
index 8a98b0a6105f62e29ded79bfe7c2a2b1a0e9f048..fafb8db5adc9adb1bfe7ce505fd58224bf3360e9 100644
--- a/tests/test_frame.py
+++ b/tests/test_frame.py
@@ -22,7 +22,7 @@ def try_remove(filename) -> None:
 # TODO add ITP parsing tests
 class FrameTest(unittest.TestCase):
     base_dir = pathlib.Path(__file__).absolute().parent
-    data_dir = base_dir.joinpath('data')
+    data_dir = base_dir.joinpath("data")
 
     def check_reference_topology(self, frame, skip_names=True):
         self.assertEqual(663, frame.natoms)
@@ -33,8 +33,8 @@ class FrameTest(unittest.TestCase):
         self.assertEqual(3, residue.n_atoms)
 
         if not skip_names:  # MDTraj renames water
-            self.assertEqual('SOL', residue.name)
-            self.assertEqual('OW', residue.atom(0).name)
+            self.assertEqual("SOL", residue.name)
+            self.assertEqual("OW", residue.atom(0).name)
 
     def check_reference_frame(self, frame, check_box: bool = True):
         self.check_reference_topology(frame)
@@ -44,122 +44,140 @@ class FrameTest(unittest.TestCase):
 
         np.testing.assert_allclose(atom0_coords, frame.atom(0).coords)
         if check_box:
-            np.testing.assert_allclose(box_vectors, frame.unitcell_lengths,
-                                       rtol=1e-4)  # PDB files are f9.3
+            np.testing.assert_allclose(
+                box_vectors, frame.unitcell_lengths, rtol=1e-4
+            )  # PDB files are f9.3
 
     def check_reference_trajectory(self, frame):
         self.check_reference_topology(frame)
-
-        atom0_coords = np.array([
-            [1.176     , 1.1520001 , 1.5860001 ],
-            [1.1220001 , 1.13      , 1.534     ],
-            [1.0580001 , 1.1620001 , 1.462     ],
-            [0.91600007, 1.276     , 1.5580001 ],
-            [0.73200005, 1.1240001 , 1.286     ],
-            [0.64000005, 1.2160001 , 1.258     ],
-            [0.632     , 1.312     , 1.2520001 ],
-            [0.606     , 1.284     , 1.246     ],
-            [0.582     , 1.312     , 1.1600001 ],
-            [0.68200004, 1.22      , 1.25      ],
-            [0.69600004, 1.33      , 1.21      ],
-        ], dtype=np.float32)  # yapf: disable
-
-        unitcell_lengths = np.array([
-            [1.9052   , 1.9052   , 1.9052   ],
-            [1.9032527, 1.9032527, 1.9032527],
-            [1.9040661, 1.9040661, 1.9040661],
-            [1.896811 , 1.896811 , 1.896811 ],
-            [1.8985983, 1.8985983, 1.8985983],
-            [1.9033976, 1.9033976, 1.9033976],
-            [1.8904614, 1.8904614, 1.8904614],
-            [1.9013108, 1.9013108, 1.9013108],
-            [1.8946321, 1.8946321, 1.8946321],
-            [1.898091 , 1.898091 , 1.898091 ],
-            [1.898684 , 1.898684 , 1.898684 ],
-        ], dtype=np.float32)  # yapf: disable
+        atom0_coords = np.array(
+            [
+                [1.176     , 1.1520001 , 1.5860001 ],
+                [1.1220001 , 1.13      , 1.534     ],
+                [1.0580001 , 1.1620001 , 1.462     ],
+                [0.91600007, 1.276     , 1.5580001 ],
+                [0.73200005, 1.1240001 , 1.286     ],
+                [0.64000005, 1.2160001 , 1.258     ],
+                [0.632     , 1.312     , 1.2520001 ],
+                [0.606     , 1.284     , 1.246     ],
+                [0.582     , 1.312     , 1.1600001 ],
+                [0.68200004, 1.22      , 1.25      ],
+                [0.69600004, 1.33      , 1.21      ],
+            ],
+            dtype=np.float32,
+        )  # fmt: skip
+        unitcell_lengths = np.array(
+            [
+                [1.9052   , 1.9052   , 1.9052   ],
+                [1.9032527, 1.9032527, 1.9032527],
+                [1.9040661, 1.9040661, 1.9040661],
+                [1.896811 , 1.896811 , 1.896811 ],
+                [1.8985983, 1.8985983, 1.8985983],
+                [1.9033976, 1.9033976, 1.9033976],
+                [1.8904614, 1.8904614, 1.8904614],
+                [1.9013108, 1.9013108, 1.9013108],
+                [1.8946321, 1.8946321, 1.8946321],
+                [1.898091 , 1.898091 , 1.898091 ],
+                [1.898684 , 1.898684 , 1.898684 ],
+            ],
+            dtype=np.float32,
+        )  # fmt: skip
 
         np.testing.assert_allclose(atom0_coords, frame.atom(0).coords)
-        np.testing.assert_allclose(unitcell_lengths,
-                                   frame.unitcell_lengths,
-                                   rtol=1e-4)  # PDB files are f9.3
+        np.testing.assert_allclose(
+            unitcell_lengths, frame.unitcell_lengths, rtol=1e-4
+        )  # PDB files are f9.3
 
     def test_frame_create(self):
         Frame()
 
     def test_frame_add_residue(self):
         frame = Frame()
-        residue = frame.add_residue('TEST_RESIDUE')
+        residue = frame.add_residue("TEST_RESIDUE")
 
         self.assertTrue(residue is frame.residue(0))
 
     def test_frame_read_gro(self):
-        frame = Frame(self.data_dir.joinpath('water.gro'))
+        frame = Frame(self.data_dir.joinpath("water.gro"))
 
         self.check_reference_frame(frame)
 
     def test_frame_read_pdb(self):
-        frame = Frame(self.data_dir.joinpath('water.pdb'))
+        frame = Frame(self.data_dir.joinpath("water.pdb"))
 
         self.check_reference_frame(frame)
 
     def test_frame_read_pdb_nobox(self):
         """Test reading a PDB with zero box volume."""
-        frame = Frame(self.data_dir.joinpath('water-nobox.pdb'))
+        frame = Frame(self.data_dir.joinpath("water-nobox.pdb"))
 
         self.check_reference_frame(frame, check_box=False)
 
     def test_frame_read_zero_box(self):
-        frame = Frame(self.data_dir.joinpath('polyethene.gro'))
+        frame = Frame(self.data_dir.joinpath("polyethene.gro"))
         self.assertIsNone(frame.unitcell_lengths)
 
     def test_frame_any_read_unsupported(self):
         with self.assertRaises(UnsupportedFormatException):
-            _ = Frame(self.data_dir.joinpath('dppc.map'))
+            _ = Frame(self.data_dir.joinpath("dppc.map"))
 
     def test_frame_output_gro(self):
-        frame = Frame(self.data_dir.joinpath('water.gro'))
+        frame = Frame(self.data_dir.joinpath("water.gro"))
 
         with tempfile.TemporaryDirectory() as tmpdir:
             tmp_path = pathlib.Path(tmpdir)
 
-            frame.save(tmp_path.joinpath('water-out.gro'))
+            frame.save(tmp_path.joinpath("water-out.gro"))
             self.assertTrue(
-                filecmp.cmp(self.data_dir.joinpath('water.gro'),
-                            tmp_path.joinpath('water-out.gro')))
+                filecmp.cmp(
+                    self.data_dir.joinpath("water.gro"),
+                    tmp_path.joinpath("water-out.gro"),
+                )
+            )
 
     def test_frame_read_xtc_numframes(self):
-        frame = Frame(self.data_dir.joinpath('water.gro'), self.data_dir.joinpath('water.xtc'))
+        frame = Frame(
+            self.data_dir.joinpath("water.gro"), self.data_dir.joinpath("water.xtc")
+        )
         self.assertEqual(11, frame.n_frames)
 
     def test_frame_read_xtc(self):
-        frame = Frame(self.data_dir.joinpath('water.gro'), self.data_dir.joinpath('water.xtc'))
+        frame = Frame(
+            self.data_dir.joinpath("water.gro"), self.data_dir.joinpath("water.xtc")
+        )
 
         self.check_reference_trajectory(frame)
 
     def test_frame_accept_path(self):
         """Test that Frame accepts :class:`pathlib.Path` arguments."""
-        _ = Frame(self.data_dir.joinpath('water.gro'),
-                  self.data_dir.joinpath('water.xtc'))
+        _ = Frame(
+            self.data_dir.joinpath("water.gro"), self.data_dir.joinpath("water.xtc")
+        )
 
     def test_frame_write_xtc(self):
         """Test that :class:`Frame` can save a trajectory file."""
-        frame = Frame(self.data_dir.joinpath('water.gro'),
-                      self.data_dir.joinpath('water.xtc'))
+        frame = Frame(
+            self.data_dir.joinpath("water.gro"), self.data_dir.joinpath("water.xtc")
+        )
 
         with tempfile.TemporaryDirectory() as tmpdir:
             tmp_path = pathlib.Path(tmpdir)
 
-            frame.save(tmp_path.joinpath('water.xtc'))
+            frame.save(tmp_path.joinpath("water.xtc"))
             self.assertTrue(
                 util.compare_trajectories(
-                    self.data_dir.joinpath('water.xtc'),
-                    tmp_path.joinpath('water.xtc'),
-                    topology_file=self.data_dir.joinpath('water.gro')))
+                    self.data_dir.joinpath("water.xtc"),
+                    tmp_path.joinpath("water.xtc"),
+                    topology_file=self.data_dir.joinpath("water.gro"),
+                )
+            )
 
     def test_raise_nonmatching_system_mdtraj(self):
         with self.assertRaises(NonMatchingSystemError):
-            _ = Frame(self.data_dir.joinpath('water.gro'), self.data_dir.joinpath('sugar.xtc'))
+            _ = Frame(
+                self.data_dir.joinpath("water.gro"), self.data_dir.joinpath("sugar.xtc")
+            )
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     unittest.main()
diff --git a/tests/test_functionalforms.py b/tests/test_functionalforms.py
index 0266386571b2ed82a8063ccb07991d95cd393b9d..f8b13d61722875f864dcf8f0e3912e08199e6434 100644
--- a/tests/test_functionalforms.py
+++ b/tests/test_functionalforms.py
@@ -29,7 +29,7 @@ class FunctionalFormTest(unittest.TestCase):
                 return "TestResultFconst"
 
         funcs = get_functional_forms()
-        test_func = funcs['TestFunc'].value
+        test_func = funcs["TestFunc"].value
 
         self.assertEqual("TestResultEqm", test_func.eqm(None, None))
         self.assertEqual("TestResultFconst", test_func.fconst(None, None))
@@ -43,18 +43,11 @@ class FunctionalFormTest(unittest.TestCase):
 
         func = funcs["Harmonic"].value()
 
-        np.testing.assert_allclose(
-            func.eqm(vals, None), np.pi
-        )
+        np.testing.assert_allclose(func.eqm(vals, None), np.pi)
 
-        func = funcs["Harmonic"].value(
-            mean_func=circular_mean
-        )
+        func = funcs["Harmonic"].value(mean_func=circular_mean)
 
-        np.testing.assert_allclose(
-            func.eqm(vals, None), 0,
-            atol=1e-6
-        )
+        np.testing.assert_allclose(func.eqm(vals, None), 0, atol=1e-6)
 
     def test_inject_variance_function(self):
         """
@@ -62,19 +55,14 @@ class FunctionalFormTest(unittest.TestCase):
         """
         funcs = get_functional_forms()
 
-        func = funcs["Harmonic"].value(
-            variance_func=circular_variance
-        )
+        func = funcs["Harmonic"].value(variance_func=circular_variance)
 
         vals = [0, 0.5 * np.pi]
         ref = func.fconst(vals, 300)
 
         vals = [0.25 * np.pi, 1.75 * np.pi]
-        np.testing.assert_allclose(
-            func.fconst(vals, 300), ref,
-            atol=0
-        )
+        np.testing.assert_allclose(func.fconst(vals, 300), ref, atol=0)
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     unittest.main()
diff --git a/tests/test_mapping.py b/tests/test_mapping.py
index b8d2d93c9c77489e8e264d197794cef1973eab5e..273f6478ce1e25db0e0c2b31c08494d0536a8256 100644
--- a/tests/test_mapping.py
+++ b/tests/test_mapping.py
@@ -17,18 +17,18 @@ class DummyOptions:
 class BeadMapTest(unittest.TestCase):
     def test_beadmap_create(self):
         """Test that a BeadMap can be created."""
-        bead = BeadMap('BEAD', 1, atoms=['ATOM1', 'ATOM2'])
+        bead = BeadMap("BEAD", 1, atoms=["ATOM1", "ATOM2"])
 
-        self.assertEqual('BEAD', bead.name)
+        self.assertEqual("BEAD", bead.name)
         self.assertEqual(1, bead.num)
 
         self.assertEqual(2, len(bead))
-        self.assertEqual('ATOM1', bead[0])
-        self.assertEqual('ATOM2', bead[1])
+        self.assertEqual("ATOM1", bead[0])
+        self.assertEqual("ATOM2", bead[1])
 
     def test_beadmap_guess_masses(self):
         """Test that masses can be assigned from atom names."""
-        bead = BeadMap('BEAD', 1, atoms=['C1', 'H1', 'H2'])
+        bead = BeadMap("BEAD", 1, atoms=["C1", "H1", "H2"])
         bead.guess_atom_masses()
 
         # Total mass of bead
@@ -36,173 +36,209 @@ class BeadMapTest(unittest.TestCase):
 
         # Contribution of each atom to bead mass
         np.testing.assert_allclose(
-            np.array([12, 1, 1]) / 14,
-            bead.weights_dict['mass'],
-            rtol=0.01)  # Account for 1% error in approximate masses
+            np.array([12, 1, 1]) / 14, bead.weights_dict["mass"], rtol=0.01
+        )  # Account for 1% error in approximate masses
 
 
 class MappingTest(unittest.TestCase):
     base_dir = pathlib.Path(__file__).absolute().parent
-    data_dir = base_dir.joinpath('data')
+    data_dir = base_dir.joinpath("data")
 
     def test_mapping_create(self):
         """Test that a mapping can be correctly read from a file."""
-        mapping = Mapping(self.data_dir.joinpath('water.map'), DummyOptions)
+        mapping = Mapping(self.data_dir.joinpath("water.map"), DummyOptions)
         self.assertEqual(2, len(mapping))  # SOL and HOH
-        self.assertTrue('SOL' in mapping)
+        self.assertTrue("SOL" in mapping)
 
-        sol_map = mapping['SOL']
+        sol_map = mapping["SOL"]
 
         self.assertEqual(1, len(sol_map))
         self.assertEqual(3, len(sol_map[0].atoms))
-        self.assertEqual('OW', sol_map[0].atoms[0])
-        self.assertEqual('HW1', sol_map[0].atoms[1])
-        self.assertEqual('HW2', sol_map[0].atoms[2])
+        self.assertEqual("OW", sol_map[0].atoms[0])
+        self.assertEqual("HW1", sol_map[0].atoms[1])
+        self.assertEqual("HW2", sol_map[0].atoms[2])
 
     def test_mapping_rename(self):
         """Test that an alternate mapping is created using MDTraj conventions."""
-        mapping = Mapping(self.data_dir.joinpath('water.map'), DummyOptions)
+        mapping = Mapping(self.data_dir.joinpath("water.map"), DummyOptions)
         self.assertEqual(2, len(mapping))  # SOL and HOH
-        self.assertTrue('HOH' in mapping)
+        self.assertTrue("HOH" in mapping)
 
-        sol_map = mapping['HOH']
+        sol_map = mapping["HOH"]
 
         self.assertEqual(1, len(sol_map))
         self.assertEqual(3, len(sol_map[0].atoms))
-        self.assertEqual('O', sol_map[0].atoms[0])
-        self.assertEqual('H1', sol_map[0].atoms[1])
-        self.assertEqual('H2', sol_map[0].atoms[2])
+        self.assertEqual("O", sol_map[0].atoms[0])
+        self.assertEqual("H1", sol_map[0].atoms[1])
+        self.assertEqual("H2", sol_map[0].atoms[2])
 
     def test_virtual_mapping_create(self):
-        mapping = Mapping(self.data_dir.joinpath('martini3/naphthalene.map'), DummyOptions)
+        mapping = Mapping(
+            self.data_dir.joinpath("martini3/naphthalene.map"), DummyOptions
+        )
         self.assertEqual(1, len(mapping))
         self.assertTrue("NAPH" in mapping)
         self.assertEqual(5, len(mapping["NAPH"]))
         self.assertTrue(isinstance(mapping["NAPH"][2], VirtualMap))
         self.assertEqual(4, len(mapping["NAPH"][2].atoms))
         self.assertEqual("R1", mapping["NAPH"][2].atoms[0])
-        self.assertEqual(1, [isinstance(bead, VirtualMap) for bead in mapping["NAPH"]].count(True))
+        self.assertEqual(
+            1, [isinstance(bead, VirtualMap) for bead in mapping["NAPH"]].count(True)
+        )
 
     def test_mapping_apply(self):
-        mapping = Mapping(self.data_dir.joinpath('water.map'), DummyOptions)
-        frame = Frame(self.data_dir.joinpath('water.gro'))
+        mapping = Mapping(self.data_dir.joinpath("water.map"), DummyOptions)
+        frame = Frame(self.data_dir.joinpath("water.gro"))
         cg_frame = mapping.apply(frame)
 
         self.assertEqual(frame.natoms / 3, cg_frame.natoms)
 
         cg_frame.save("water-cg.gro")
 
-        self.assertTrue(filecmp.cmp(self.data_dir.joinpath('water-cg.gro'), "water-cg.gro"))
+        self.assertTrue(
+            filecmp.cmp(self.data_dir.joinpath("water-cg.gro"), "water-cg.gro")
+        )
         os.remove("water-cg.gro")
 
     def test_mapping_duplicate_atom_weight(self):
         """Test that an atom can be duplicated in a bead specification to act as a weighting."""
-        frame = Frame(self.data_dir.joinpath('water.gro'))
+        frame = Frame(self.data_dir.joinpath("water.gro"))
 
-        mapping = Mapping(self.data_dir.joinpath('water.map'), DummyOptions)
+        mapping = Mapping(self.data_dir.joinpath("water.map"), DummyOptions)
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([[0.716, 1.300, 1.198]]),
-                                   cg_frame.residue(0).atom(0).coords, atol=0.0005)
-
-        mapping = Mapping(self.data_dir.joinpath('water_double_weight.map'), DummyOptions)
+        np.testing.assert_allclose(
+            np.array([[0.716, 1.300, 1.198]]),
+            cg_frame.residue(0).atom(0).coords,
+            atol=0.0005,
+        )
+
+        mapping = Mapping(
+            self.data_dir.joinpath("water_double_weight.map"), DummyOptions
+        )
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([[0.720, 1.294, 1.195]]),
-                                   cg_frame.residue(0).atom(0).coords, atol=0.0005)
+        np.testing.assert_allclose(
+            np.array([[0.720, 1.294, 1.195]]),
+            cg_frame.residue(0).atom(0).coords,
+            atol=0.0005,
+        )
 
     def test_mapping_charges(self):
-        mapping = Mapping(self.data_dir.joinpath('dppc.map'), DummyOptions)
+        mapping = Mapping(self.data_dir.joinpath("dppc.map"), DummyOptions)
         self.assertEqual(1, mapping["DPPC"][0].charge)
         self.assertEqual(-1, mapping["DPPC"][1].charge)
 
     def test_mapping_pbc(self):
-        frame = Frame(self.data_dir.joinpath('pbcwater.gro'))
+        frame = Frame(self.data_dir.joinpath("pbcwater.gro"))
 
-        mapping = Mapping(self.data_dir.joinpath('water.map'), DummyOptions)
+        mapping = Mapping(self.data_dir.joinpath("water.map"), DummyOptions)
         cg_frame = mapping.apply(frame)
 
         np.testing.assert_allclose(frame.atom(0).coords, cg_frame.atom(0).coords)
 
     def test_mapping_weights_geom(self):
-        frame = Frame(self.data_dir.joinpath('two.gro'))
+        frame = Frame(self.data_dir.joinpath("two.gro"))
 
-        mapping = Mapping(self.data_dir.joinpath('two.map'), DummyOptions)
+        mapping = Mapping(self.data_dir.joinpath("two.map"), DummyOptions)
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([[1.5, 1.5, 1.5]]),
-                                   cg_frame.residue(0).atom(0).coords)
+        np.testing.assert_allclose(
+            np.array([[1.5, 1.5, 1.5]]), cg_frame.residue(0).atom(0).coords
+        )
 
     def test_virtual_mapping_weights_geom(self):
-        frame = Frame(self.data_dir.joinpath('martini3/four.gro'))
+        frame = Frame(self.data_dir.joinpath("martini3/four.gro"))
 
-        mapping = Mapping(self.data_dir.joinpath('martini3/four.map'), DummyOptions)
+        mapping = Mapping(self.data_dir.joinpath("martini3/four.map"), DummyOptions)
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([[2.5, 2.5, 2.5]]),
-                                   cg_frame.residue(0).atom(2).coords)
+        np.testing.assert_allclose(
+            np.array([[2.5, 2.5, 2.5]]), cg_frame.residue(0).atom(2).coords
+        )
 
     def test_mapping_weights_mass(self):
-        frame = Frame(self.data_dir.joinpath('two.gro'))
+        frame = Frame(self.data_dir.joinpath("two.gro"))
         options = DummyOptions()
         options.map_center = "mass"
 
-        mapping = Mapping(self.data_dir.joinpath('two.map'), options, itp_filename=self.data_dir.joinpath('two.itp'))
+        mapping = Mapping(
+            self.data_dir.joinpath("two.map"),
+            options,
+            itp_filename=self.data_dir.joinpath("two.itp"),
+        )
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([[2., 2., 2.]]),
-                                   cg_frame.residue(0).atom(0).coords)
+        np.testing.assert_allclose(
+            np.array([[2.0, 2.0, 2.0]]), cg_frame.residue(0).atom(0).coords
+        )
 
     def test_virtual_mapping_weights_mass(self):
-        frame = Frame(self.data_dir.joinpath('martini3/four.gro'))
+        frame = Frame(self.data_dir.joinpath("martini3/four.gro"))
         options = DummyOptions()
         options.virtual_map_center = "mass"
 
-        mapping = Mapping(self.data_dir.joinpath('martini3/four.map'),
-                          options,
-                          itp_filename=self.data_dir.joinpath('martini3/four.itp'))
+        mapping = Mapping(
+            self.data_dir.joinpath("martini3/four.map"),
+            options,
+            itp_filename=self.data_dir.joinpath("martini3/four.itp"),
+        )
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([[3.0, 3.0, 3.0]]),
-                                   cg_frame.residue(0).atom(2).coords)
+        np.testing.assert_allclose(
+            np.array([[3.0, 3.0, 3.0]]), cg_frame.residue(0).atom(2).coords
+        )
 
     def test_mapping_weights_guess_mass(self):
-        frame = Frame(self.data_dir.joinpath('two.gro'))
+        frame = Frame(self.data_dir.joinpath("two.gro"))
         options = DummyOptions()
         options.map_center = "mass"
 
-        mapping = Mapping(self.data_dir.joinpath('two.map'), options)
+        mapping = Mapping(self.data_dir.joinpath("two.map"), options)
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([[1.922575,  1.922575,  1.922575]], dtype=np.float32),
-                                   cg_frame.residue(0).atom(0).coords, rtol=0.01)
+        np.testing.assert_allclose(
+            np.array([[1.922575, 1.922575, 1.922575]], dtype=np.float32),
+            cg_frame.residue(0).atom(0).coords,
+            rtol=0.01,
+        )
 
     def test_virtual_mapping_weights_guess_mass(self):
-        frame = Frame(self.data_dir.joinpath('martini3/four.gro'))
+        frame = Frame(self.data_dir.joinpath("martini3/four.gro"))
         options = DummyOptions()
         options.virtual_map_center = "mass"
 
-        mapping = Mapping(self.data_dir.joinpath('martini3/four.map'), options)
+        mapping = Mapping(self.data_dir.joinpath("martini3/four.map"), options)
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([[2.83337, 2.83337, 2.83337]], dtype=np.float32),
-                                   cg_frame.residue(0).atom(2).coords, rtol=0.01)
+        np.testing.assert_allclose(
+            np.array([[2.83337, 2.83337, 2.83337]], dtype=np.float32),
+            cg_frame.residue(0).atom(2).coords,
+            rtol=0.01,
+        )
 
     def test_mapping_weights_first(self):
-        frame = Frame(self.data_dir.joinpath('two.gro'))
+        frame = Frame(self.data_dir.joinpath("two.gro"))
         options = DummyOptions()
         options.map_center = "first"
 
-        mapping = Mapping(self.data_dir.joinpath('two.map'), options, itp_filename=self.data_dir.joinpath('two.itp'))
+        mapping = Mapping(
+            self.data_dir.joinpath("two.map"),
+            options,
+            itp_filename=self.data_dir.joinpath("two.itp"),
+        )
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([[1., 1., 1.]]),
-                                   cg_frame.residue(0).atom(0).coords)
+        np.testing.assert_allclose(
+            np.array([[1.0, 1.0, 1.0]]), cg_frame.residue(0).atom(0).coords
+        )
 
     def test_mapping_itp_multi(self):
-        mapping = Mapping(self.data_dir.joinpath('membrane/membrane.map'),
-                          DummyOptions,
-                          itp_filename=self.data_dir.joinpath('membrane/membrane.top'))
+        mapping = Mapping(
+            self.data_dir.joinpath("membrane/membrane.map"),
+            DummyOptions,
+            itp_filename=self.data_dir.joinpath("membrane/membrane.top"),
+        )
         self.assertAlmostEqual(-1.2, mapping["POPE"][0].charge, delta=0.0001)
         self.assertAlmostEqual(0, mapping["POPG"][0].charge, delta=0.0001)
 
@@ -213,8 +249,8 @@ class MappingTest(unittest.TestCase):
 
         This is instead checked and logged as a warning in the outer script.
         """
-        frame = Frame(self.data_dir.joinpath('sugar.gro'))
-        mapping = Mapping(self.data_dir.joinpath('water.map'), DummyOptions)
+        frame = Frame(self.data_dir.joinpath("sugar.gro"))
+        mapping = Mapping(self.data_dir.joinpath("water.map"), DummyOptions)
 
         cg_frame = mapping.apply(frame)
 
diff --git a/tests/test_parsers_cfg.py b/tests/test_parsers_cfg.py
index 68e7760b7df5daa118bec4cced856c36eabe42ec..ea20a4a9a2b7a0a044979bdaa4c88078a6d51743 100644
--- a/tests/test_parsers_cfg.py
+++ b/tests/test_parsers_cfg.py
@@ -7,46 +7,46 @@ from pycgtool.parsers import CFG
 
 class TestParsersCFG(unittest.TestCase):
     base_dir = pathlib.Path(__file__).absolute().parent
-    data_dir = base_dir.joinpath('data')
+    data_dir = base_dir.joinpath("data")
 
-    water_name = 'SOL'
+    water_name = "SOL"
     watline = ("W", "P4", "OW", "HW1", "HW2")
 
     def test_cfg_with(self):
-        with CFG(self.data_dir.joinpath('water.map')):
+        with CFG(self.data_dir.joinpath("water.map")):
             pass
 
     def test_cfg_iterate_sections(self):
-        with CFG(self.data_dir.joinpath('water.map')) as cfg:
+        with CFG(self.data_dir.joinpath("water.map")) as cfg:
             for name in cfg:
                 self.assertEqual(self.water_name, name)
 
     def test_cfg_iterate_lines(self):
-        with CFG(self.data_dir.joinpath('water.map')) as cfg:
+        with CFG(self.data_dir.joinpath("water.map")) as cfg:
             for name, section in cfg.items():
                 self.assertEqual(self.water_name, name)
                 for line in section:
                     self.assertEqual(self.watline, line)
 
     def test_cfg_get_section(self):
-        with CFG(self.data_dir.joinpath('water.map')) as cfg:
+        with CFG(self.data_dir.joinpath("water.map")) as cfg:
             self.assertTrue(self.water_name in cfg)
             for line in cfg[self.water_name]:
                 self.assertEqual(self.watline, line)
 
     def test_cfg_duplicate_error(self):
         with self.assertRaises(DuplicateSectionError):
-            CFG(self.data_dir.joinpath('twice.cfg'))
+            CFG(self.data_dir.joinpath("twice.cfg"))
 
     def test_include_file(self):
-        with CFG(self.data_dir.joinpath('martini.map')) as cfg:
+        with CFG(self.data_dir.joinpath("martini.map")) as cfg:
             self.assertTrue("DOPC" in cfg)
             self.assertTrue("GLY" in cfg)
 
     def test_error_no_sections(self):
         with self.assertRaises(NoSectionError):
-            CFG(self.data_dir.joinpath('nosections.cfg'))
+            CFG(self.data_dir.joinpath("nosections.cfg"))
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     unittest.main()
diff --git a/tests/test_parsers_topology.py b/tests/test_parsers_topology.py
index 99133f293bc38d865f7a695bf37bf9a6e2e42324..289230d952baa73e311e915000ff22e9bf125f18 100644
--- a/tests/test_parsers_topology.py
+++ b/tests/test_parsers_topology.py
@@ -7,40 +7,40 @@ from pycgtool.parsers import ITP
 
 class TestParsersITP(unittest.TestCase):
     base_dir = pathlib.Path(__file__).absolute().parent
-    data_dir = base_dir.joinpath('data')
+    data_dir = base_dir.joinpath("data")
 
     def test_empty_itp(self):
         itp = ITP()
         self.assertEqual(0, len(itp))
 
     def test_itp_with(self):
-        with ITP(self.data_dir.joinpath('membrane/membrane.top')):
+        with ITP(self.data_dir.joinpath("membrane/membrane.top")):
             pass
 
     def test_itp_duplicate_error(self):
         with self.assertRaises(DuplicateSectionError):
-            ITP(self.data_dir.joinpath('membrane/twice.itp'))
+            ITP(self.data_dir.joinpath("membrane/twice.itp"))
 
     def test_itp_include(self):
-        with ITP(self.data_dir.joinpath('membrane/membrane.top')) as itp:
-            self.assertIn('CDL2', itp)
-            self.assertIn('POPE', itp)
+        with ITP(self.data_dir.joinpath("membrane/membrane.top")) as itp:
+            self.assertIn("CDL2", itp)
+            self.assertIn("POPE", itp)
 
     def test_read_section(self):
         pii_line = ["1", "CC3161", "1", "PII", "C13", "1", "0.1400", "12.0110"]
 
-        with ITP(self.data_dir.joinpath('membrane/membrane.top')) as itp:
+        with ITP(self.data_dir.joinpath("membrane/membrane.top")) as itp:
             self.assertEqual(len(itp["PII"]["atoms"]), 283)
             self.assertListEqual(list(itp["PII"]["atoms"][0]), pii_line)
 
     def test_double_sections(self):
-        with ITP(self.data_dir.joinpath('two_double_section.itp')) as itp:
+        with ITP(self.data_dir.joinpath("two_double_section.itp")) as itp:
             self.assertEqual(len(itp["TWO"]["atoms"]), 2)
 
     def test_unknown_sections(self):
-        with ITP(self.data_dir.joinpath('two_unknown_section.itp')) as itp:
+        with ITP(self.data_dir.joinpath("two_unknown_section.itp")) as itp:
             self.assertEqual(len(itp["TWO"]["atoms"]), 2)
 
     def test_error_no_sections(self):
         with self.assertRaises(NoSectionError):
-            ITP(self.data_dir.joinpath('nosections.cfg'))
+            ITP(self.data_dir.joinpath("nosections.cfg"))
diff --git a/tests/test_pycgtool.py b/tests/test_pycgtool.py
index bfecb005936f3cb0d1dec1dc36c32dc2c2788a6c..8eb365cbcf010422bbce60c5fb5dc87f62e7f635 100644
--- a/tests/test_pycgtool.py
+++ b/tests/test_pycgtool.py
@@ -13,22 +13,24 @@ import pycgtool.__main__ as main
 
 def get_args(name, out_dir, extra: typing.Optional[typing.Mapping] = None):
     base_dir = pathlib.Path(__file__).absolute().parent
-    data_dir = base_dir.joinpath('data')
+    data_dir = base_dir.joinpath("data")
 
     args = [
-        data_dir.joinpath(f'{name}.gro'),
-        data_dir.joinpath(f'{name}.xtc'),
+        data_dir.joinpath(f"{name}.gro"),
+        data_dir.joinpath(f"{name}.xtc"),
     ]
 
     kwargs = {
-        '-m': data_dir.joinpath(f'{name}.map'),
-        '-b': data_dir.joinpath(f'{name}.bnd'),
-        '--out-dir': out_dir,
+        "-m": data_dir.joinpath(f"{name}.map"),
+        "-b": data_dir.joinpath(f"{name}.bnd"),
+        "--out-dir": out_dir,
     }
 
     parsed_args = main.parse_arguments(
-        itertools.chain(map(str, args),
-                        *[[key, str(value)] for key, value in kwargs.items()]))
+        itertools.chain(
+            map(str, args), *[[key, str(value)] for key, value in kwargs.items()]
+        )
+    )
 
     if extra is not None:
         for key, value in extra.items():
@@ -40,36 +42,42 @@ def get_args(name, out_dir, extra: typing.Optional[typing.Mapping] = None):
 
 class PycgtoolTest(unittest.TestCase):
     base_dir = pathlib.Path(__file__).absolute().parent
-    data_dir = base_dir.joinpath('data')
+    data_dir = base_dir.joinpath("data")
 
     def test_pycgtool_help(self):
         subprocess.check_call(["python", "-m", "pycgtool", "-h"])
         subprocess.check_call(["python", "-m", "pycgtool", "--help"])
 
     def test_pycgtool_version(self):
-        subprocess.check_call(['python', '-m', 'pycgtool', '-v'])
-        subprocess.check_call(['python', '-m', 'pycgtool', '--version'])
+        subprocess.check_call(["python", "-m", "pycgtool", "-v"])
+        subprocess.check_call(["python", "-m", "pycgtool", "--version"])
 
     def test_parse_arguments(self):
-        args = main.parse_arguments([
-            'TOPOLOGY',
-            '-m',
-            'MAP',
-            '--begin',
-            '1000',
-        ])
-
-        self.assertEqual('TOPOLOGY', args.topology)
-        self.assertEqual('MAP', args.mapping)
+        args = main.parse_arguments(
+            [
+                "TOPOLOGY",
+                "-m",
+                "MAP",
+                "--begin",
+                "1000",
+            ]
+        )
+
+        self.assertEqual("TOPOLOGY", args.topology)
+        self.assertEqual("MAP", args.mapping)
         self.assertEqual(1000, args.begin)
 
     def test_map_only(self):
         with tempfile.TemporaryDirectory() as tmpdir:
             tmp_path = pathlib.Path(tmpdir)
-            args = get_args('sugar', tmp_path, extra={
-                'output_xtc': True,
-                'bnd': None,
-            })
+            args = get_args(
+                "sugar",
+                tmp_path,
+                extra={
+                    "output_xtc": True,
+                    "bnd": None,
+                },
+            )
 
             # Equivalent to
             # pycgtool <top> <trj> -m <map> --output-xtc
@@ -77,14 +85,16 @@ class PycgtoolTest(unittest.TestCase):
 
             self.assertTrue(
                 util.compare_trajectories(
-                    self.data_dir.joinpath('sugar_out.xtc'),
-                    tmp_path.joinpath('out.xtc'),
-                    topology_file=self.data_dir.joinpath('sugar_out.gro')))
+                    self.data_dir.joinpath("sugar_out.xtc"),
+                    tmp_path.joinpath("out.xtc"),
+                    topology_file=self.data_dir.joinpath("sugar_out.gro"),
+                )
+            )
 
     def test_full(self):
         with tempfile.TemporaryDirectory() as tmpdir:
             tmp_path = pathlib.Path(tmpdir)
-            args = get_args('sugar', tmp_path)
+            args = get_args("sugar", tmp_path)
 
             # Equivalent to
             # pycgtool <top> <trj> -m <map> -b <bnd>
@@ -95,87 +105,103 @@ class PycgtoolTest(unittest.TestCase):
                     tmp_path.joinpath("out.itp"),
                     self.data_dir.joinpath("sugar_out.itp"),
                     rtol=0.001,
-                    verbose=True))
+                    verbose=True,
+                )
+            )
 
             self.assertTrue(
                 util.compare_trajectories(
-                    self.data_dir.joinpath('sugar_out.gro'),
-                    tmp_path.joinpath('out.gro'),
-                    rtol=0.001))
+                    self.data_dir.joinpath("sugar_out.gro"),
+                    tmp_path.joinpath("out.gro"),
+                    rtol=0.001,
+                )
+            )
 
     def test_full_no_traj(self):
         with tempfile.TemporaryDirectory() as tmpdir:
             tmp_path = pathlib.Path(tmpdir)
-            args = get_args('sugar', tmp_path, extra={
-                'trajectory': None,
-            })
+            args = get_args(
+                "sugar",
+                tmp_path,
+                extra={
+                    "trajectory": None,
+                },
+            )
 
             # Equivalent to
             # pycgtool <top> -m <map> -b <bnd>
             main.PyCGTOOL(args)
 
-            self.assertFalse(tmp_path.joinpath('out.itp').exists())
+            self.assertFalse(tmp_path.joinpath("out.itp").exists())
 
             self.assertTrue(
                 util.compare_trajectories(
-                    self.data_dir.joinpath('sugar_out.gro'),
-                    tmp_path.joinpath('out.gro'),
-                    rtol=0.001))
+                    self.data_dir.joinpath("sugar_out.gro"),
+                    tmp_path.joinpath("out.gro"),
+                    rtol=0.001,
+                )
+            )
 
     def test_measure_only(self):
         with tempfile.TemporaryDirectory() as tmpdir:
             tmp_path = pathlib.Path(tmpdir)
-            args = get_args('sugar', tmp_path, extra={
-                'mapping': None,
-                'trajectory': None,
-            })
+            args = get_args(
+                "sugar",
+                tmp_path,
+                extra={
+                    "mapping": None,
+                    "trajectory": None,
+                },
+            )
 
             # Equivalent to
             # pycgtool <top> -b <bnd>
             main.PyCGTOOL(args)
 
             # Does not produce itp file
-            self.assertFalse(tmp_path.joinpath('out.itp').exists())
+            self.assertFalse(tmp_path.joinpath("out.itp").exists())
 
             # Check bond dump files against reference
-            for bond_type in ['length', 'angle', 'dihedral']:
-                out_file = tmp_path.joinpath(f'ALLA_{bond_type}.dat')
+            for bond_type in ["length", "angle", "dihedral"]:
+                out_file = tmp_path.joinpath(f"ALLA_{bond_type}.dat")
                 self.assertTrue(out_file.exists())
 
-                self.assertTrue(util.cmp_file_whitespace_float(
-                    self.data_dir.joinpath(f'ALLA_{bond_type}_one.dat'),
-                    out_file
-                ))
+                self.assertTrue(
+                    util.cmp_file_whitespace_float(
+                        self.data_dir.joinpath(f"ALLA_{bond_type}_one.dat"), out_file
+                    )
+                )
 
     def test_forcefield(self):
         with tempfile.TemporaryDirectory() as tmpdir:
             tmp_path = pathlib.Path(tmpdir)
-            args = get_args('sugar',
-                            tmp_path,
-                            extra={
-                                'output_forcefield': True,
-                            })
+            args = get_args(
+                "sugar",
+                tmp_path,
+                extra={
+                    "output_forcefield": True,
+                },
+            )
 
             # Equivalent to
             # pycgtool <top> <trj> -m <map> -b <bnd> --output-forcefield
             main.PyCGTOOL(args)
 
             # Does not produce itp file
-            self.assertFalse(tmp_path.joinpath('out.itp').exists())
+            self.assertFalse(tmp_path.joinpath("out.itp").exists())
 
-            out_ff_dir = tmp_path.joinpath('ffout.ff')
+            out_ff_dir = tmp_path.joinpath("ffout.ff")
             self.assertTrue(out_ff_dir.is_dir())
 
             # Compare all files in ffout.ff to reference versions
             for out_file in out_ff_dir.iterdir():
                 ref_file = self.data_dir.joinpath(out_file)
 
-                self.assertTrue(
-                    util.cmp_file_whitespace_float(ref_file, out_file))
+                self.assertTrue(util.cmp_file_whitespace_float(ref_file, out_file))
 
     # TODO more tests
     # TODO test wrong args.end
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     unittest.main()
diff --git a/tests/test_util.py b/tests/test_util.py
index 76fc44b1d70a114a63254b2ed38887a0995230a6..2015c3a77394dfb6b479be027d0de464f76b3b86 100644
--- a/tests/test_util.py
+++ b/tests/test_util.py
@@ -49,9 +49,9 @@ class UtilTest(unittest.TestCase):
             tmp_dir = pathlib.Path(t)
 
             paths = [
-                tmp_dir.joinpath('testfile'),
-                tmp_dir.joinpath('#testfile.1#'),
-                tmp_dir.joinpath('#testfile.2#'),
+                tmp_dir.joinpath("testfile"),
+                tmp_dir.joinpath("#testfile.1#"),
+                tmp_dir.joinpath("#testfile.2#"),
             ]
             path = paths[0]
 
@@ -85,9 +85,9 @@ class UtilTest(unittest.TestCase):
             tmp_dir = pathlib.Path(t)
 
             paths = [
-                tmp_dir.joinpath('testfile'),
-                tmp_dir.joinpath('#testfile.1#'),
-                tmp_dir.joinpath('#testfile.2#'),
+                tmp_dir.joinpath("testfile"),
+                tmp_dir.joinpath("#testfile.1#"),
+                tmp_dir.joinpath("#testfile.2#"),
             ]
             path = paths[0]
 
@@ -168,27 +168,31 @@ class UtilTest(unittest.TestCase):
 
     def test_cmp_whitespace_float_diff_lengths(self):
         """Test that lists of different lengths are considered different."""
-        ref = ['1']
-        test = ['1', '2']
+        ref = ["1"]
+        test = ["1", "2"]
 
         self.assertFalse(cmp_whitespace_float(ref, test))
 
     def test_circular_mean(self):
         values = np.deg2rad(np.array([-175, 165], dtype=np.float32))
-        test = np.deg2rad(175.)
-        self.assertAlmostEqual(circular_mean(values),  test, delta=test / 500.)
+        test = np.deg2rad(175.0)
+        self.assertAlmostEqual(circular_mean(values), test, delta=test / 500.0)
 
         values = np.deg2rad(np.array([165, 145], dtype=np.float32))
-        test = np.deg2rad(155.)
-        self.assertAlmostEqual(circular_mean(values), test, delta=test / 500.)
+        test = np.deg2rad(155.0)
+        self.assertAlmostEqual(circular_mean(values), test, delta=test / 500.0)
 
     def test_circular_variance(self):
         values = np.deg2rad(np.array([165, 145], dtype=np.float32))
         test_var = np.var(values)
-        self.assertAlmostEqual(circular_variance(values), test_var, delta=test_var / 500.)
+        self.assertAlmostEqual(
+            circular_variance(values), test_var, delta=test_var / 500.0
+        )
 
         values = np.deg2rad(np.array([-175, 165], dtype=np.float32))
-        self.assertAlmostEqual(circular_variance(values), test_var, delta=test_var / 500.)
+        self.assertAlmostEqual(
+            circular_variance(values), test_var, delta=test_var / 500.0
+        )
 
 
 class UtilFileWriteLinesTest(unittest.TestCase):
@@ -224,36 +228,38 @@ class UtilFileWriteLinesTest(unittest.TestCase):
 
 class CompareTrajectoryTest(unittest.TestCase):
     base_dir = pathlib.Path(__file__).absolute().parent
-    data_dir = base_dir.joinpath('data')
+    data_dir = base_dir.joinpath("data")
 
     def test_compare_trajectory_single(self):
-        self.assertTrue(util.compare_trajectories(
-            self.data_dir.joinpath('sugar.gro'),
-            self.data_dir.joinpath('sugar.gro')
-        ))
+        self.assertTrue(
+            util.compare_trajectories(
+                self.data_dir.joinpath("sugar.gro"), self.data_dir.joinpath("sugar.gro")
+            )
+        )
 
     def test_compare_trajectory_single_false(self):
         with self.assertRaises(ValueError):
             util.compare_trajectories(
-                self.data_dir.joinpath('sugar.gro'),
-                self.data_dir.joinpath('water.gro')
+                self.data_dir.joinpath("sugar.gro"), self.data_dir.joinpath("water.gro")
             )
 
     def test_compare_trajectory(self):
-        self.assertTrue(util.compare_trajectories(
-            self.data_dir.joinpath('sugar.xtc'),
-            self.data_dir.joinpath('sugar.xtc'),
-            topology_file=self.data_dir.joinpath('sugar.gro')
-        ))
+        self.assertTrue(
+            util.compare_trajectories(
+                self.data_dir.joinpath("sugar.xtc"),
+                self.data_dir.joinpath("sugar.xtc"),
+                topology_file=self.data_dir.joinpath("sugar.gro"),
+            )
+        )
 
     def test_compare_trajectory_false(self):
         with self.assertRaises(ValueError):
             util.compare_trajectories(
-                self.data_dir.joinpath('sugar.xtc'),
-                self.data_dir.joinpath('water.xtc'),
-                topology_file=self.data_dir.joinpath('sugar.gro')
+                self.data_dir.joinpath("sugar.xtc"),
+                self.data_dir.joinpath("water.xtc"),
+                topology_file=self.data_dir.joinpath("sugar.gro"),
             )
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     unittest.main()
diff --git a/tox.ini b/tox.ini
index ee519819c0d54f96e76e6c839f382e464edb6614..6a14acffd70fad11f629fa6bf4477da950ec9c34 100644
--- a/tox.ini
+++ b/tox.ini
@@ -4,8 +4,10 @@ envlist = lint,py3
 
 [testenv:lint]
 deps =
+    black
     flake8
 commands =
+    black --check pycgtool/ tests/
     flake8 pycgtool/ tests/
 
 [testenv]