diff --git a/create_bib_file.py b/create_bib_file.py
index 607cc3dddf7c30faa2e26fa43a7e17ebda5353ee..f70eb15c8c0b3e4be5c91d9f17dde3bed904044c 100644
--- a/create_bib_file.py
+++ b/create_bib_file.py
@@ -9,7 +9,7 @@ import requests
 import yaml
 
 
-def replace_key(key, bib_entry):
+def edit_raw_bibtex_entry(key, bib_entry):
     bib_type, *_ = bib_entry.split("{")
     _, *rest = bib_entry.split(",")
     rest = ",".join(rest)
@@ -40,8 +40,6 @@ def replace_key(key, bib_entry):
         ),  # fix for PhysRevLett.96.026804
     ]
 
-    # I got these by using JabRef and converting to abbr journals
-    # and parsing the git diff.
     journals = [
         ("Advanced Materials", "Adv. Mater."),
         ("Annals of Physics", "Ann. Phys."),
@@ -73,10 +71,6 @@ def replace_key(key, bib_entry):
         ("Science Advances", "Sci. Adv."),
         ("Scientific Reports", "Sci. Rep."),
         ("Semiconductor Science and Technology", "Semicond. Sci. Technol."),
-    ]
-
-    # Manually added
-    journals += [
         (
             "Annual Review of Condensed Matter Physics",
             "Annu. Rev. Condens. Matter Phys.",
@@ -108,12 +102,12 @@ def doi2bib(doi):
     return r.text
 
 
-fname = 'paper.yaml'
+fname = "paper.yaml"
 print("Reading: ", fname)
 
 with open(fname) as f:
-    mapping = yaml.safe_load(f)
-dois = dict(sorted(mapping.items()))
+    dois = yaml.safe_load(f)
+dois = dict(sorted(dois.items()))
 
 
 with ThreadPoolExecutor() as ex:
@@ -121,19 +115,19 @@ with ThreadPoolExecutor() as ex:
     bibs = list(futs)
 
 
-entries = [replace_key(key, bib) for key, bib in zip(dois.keys(), bibs)]
+entries = [edit_raw_bibtex_entry(key, bib) for key, bib in zip(dois.keys(), bibs)]
 
 
-with open("paper.bib", "w") as outfile:
+with open("paper.bib", "w") as out_file:
     fname = "not_on_crossref.bib"
-    outfile.write("@preamble{ {\\providecommand{\\BIBYu}{Yu} } }\n\n")
-    outfile.write(f"\n% Below is from `{fname}`.\n\n")
-    with open(fname) as infile:
-        outfile.write(infile.read())
-    outfile.write("\n% Below is from `paper.yaml` files.\n\n")
+    out_file.write("@preamble{ {\\providecommand{\\BIBYu}{Yu} } }\n\n")
+    out_file.write(f"\n% Below is from `{fname}`.\n\n")
+    with open(fname) as in_file:
+        out_file.write(in_file.read())
+    out_file.write("\n% Below is from `paper.yaml`.\n\n")
     for e in entries:
-        for line in e.split('\n'):
+        for line in e.split("\n"):
             # Remove the url line
             if "url = {" not in line:
-                outfile.write(f"{line}\n")
-        outfile.write('\n')
+                out_file.write(f"{line}\n")
+        out_file.write("\n")
diff --git a/ipynb_filter.py b/ipynb_filter.py
index 7165a27c72d2e207412be7e6b83585bc2fc69cf7..4d8c9fbd7f299add031cb57cfde3f059cf8dd4b0 100644
--- a/ipynb_filter.py
+++ b/ipynb_filter.py
@@ -16,14 +16,16 @@ from nbconvert.preprocessors import Preprocessor
 
 class RemoveMetadata(Preprocessor):
     def preprocess(self, nb, resources):
-        nb.metadata = {"language_info": {"name":"python",
-                                         "pygments_lexer": "ipython3"}}
+        nb.metadata = {
+            "language_info": {"name": "python", "pygments_lexer": "ipython3"}
+        }
         return nb, resources
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     # The filter is getting activated
     import os
+
     git_cmd = 'git config filter.ipynb_filter.clean "jupyter nbconvert --to notebook --config ipynb_filter.py --stdin --stdout"'
     os.system(git_cmd)
 else:
@@ -31,4 +33,8 @@ else:
     c.Exporter.preprocessors = [RemoveMetadata]
     c.ClearOutputPreprocessor.enabled = True
     c.ClearOutputPreprocessor.remove_metadata_fields = [
-        "deletable", "editable", "collapsed", "scrolled"]
+        "deletable",
+        "editable",
+        "collapsed",
+        "scrolled",
+    ]
diff --git a/pfaffian.py b/pfaffian.py
index 88d3f849a5d5d196285cbe7619eff1bb0e1e84cb..651dc4e2b366f7a98fcafd77e45122992bd4ee7f 100644
--- a/pfaffian.py
+++ b/pfaffian.py
@@ -1,11 +1,12 @@
 """A package for computing Pfaffians"""
 
 
+import cmath
+import math
+
 import numpy as np
 import scipy.linalg as la
 import scipy.sparse as sp
-import math
-import cmath
 
 
 def householder_real(x):
@@ -24,7 +25,7 @@ def householder_real(x):
     if sigma == 0:
         return (np.zeros(x.shape[0]), 0, x[0])
     else:
-        norm_x = math.sqrt(x[0]**2+sigma)
+        norm_x = math.sqrt(x[0] ** 2 + sigma)
 
         v = x.copy()
 
@@ -57,17 +58,17 @@ def householder_complex(x):
     if sigma == 0:
         return (np.zeros(x.shape[0]), 0, x[0])
     else:
-        norm_x = cmath.sqrt(x[0].conjugate()*x[0]+sigma)
+        norm_x = cmath.sqrt(x[0].conjugate() * x[0] + sigma)
 
         v = x.copy()
 
-        phase = cmath.exp(1j*math.atan2(x[0].imag, x[0].real))
+        phase = cmath.exp(1j * math.atan2(x[0].imag, x[0].real))
 
-        v[0] += phase*norm_x
+        v[0] += phase * norm_x
 
         v /= np.linalg.norm(v)
 
-    return (v, 2, -phase*norm_x)
+    return (v, 2, -phase * norm_x)
 
 
 def skew_tridiagonalize(A, overwrite_a=False, calc_q=True):
@@ -90,7 +91,7 @@ def skew_tridiagonalize(A, overwrite_a=False, calc_q=True):
     # Check if matrix is square
     assert A.shape[0] == A.shape[1] > 0
     # Check if it's skew-symmetric
-    assert abs((A+A.T).max()) < 1e-14
+    assert abs((A + A.T).max()) < 1e-14
 
     n = A.shape[0]
     A = np.asarray(A)  # the slice views work only properly for arrays
@@ -109,24 +110,24 @@ def skew_tridiagonalize(A, overwrite_a=False, calc_q=True):
     if calc_q:
         Q = np.eye(A.shape[0], dtype=A.dtype)
 
-    for i in range(A.shape[0]-2):
+    for i in range(A.shape[0] - 2):
         # Find a Householder vector to eliminate the i-th column
-        v, tau, alpha = householder(A[i+1:, i])
-        A[i+1, i] = alpha
-        A[i, i+1] = -alpha
-        A[i+2:, i] = 0
-        A[i, i+2:] = 0
+        v, tau, alpha = householder(A[i + 1 :, i])
+        A[i + 1, i] = alpha
+        A[i, i + 1] = -alpha
+        A[i + 2 :, i] = 0
+        A[i, i + 2 :] = 0
 
         # Update the matrix block A(i+1:N,i+1:N)
-        w = tau * A[i+1:, i+1:] @ v.conj()
-        A[i+1:, i+1:] += np.outer(v, w)-np.outer(w, v)
+        w = tau * A[i + 1 :, i + 1 :] @ v.conj()
+        A[i + 1 :, i + 1 :] += np.outer(v, w) - np.outer(w, v)
 
         if calc_q:
             # Accumulate the individual Householder reflections
             # Accumulate them in the form P_1*P_2*..., which is
             # (..*P_2*P_1)^dagger
-            y = tau * Q[:, i+1:] @ v
-            Q[:, i+1:] -= np.outer(y, v.conj())
+            y = tau * Q[:, i + 1 :] @ v
+            Q[:, i + 1 :] -= np.outer(y, v.conj())
 
     if calc_q:
         return (np.asmatrix(A), np.asmatrix(Q))
@@ -150,7 +151,7 @@ def skew_LTL(A, overwrite_a=False, calc_L=True, calc_P=True):
     # Check if matrix is square
     assert A.shape[0] == A.shape[1] > 0
     # Check if it's skew-symmetric
-    assert abs((A+A.T).max()) < 1e-14
+    assert abs((A + A.T).max()) < 1e-14
 
     n = A.shape[0]
     A = np.asarray(A)  # the slice views work only properly for arrays
@@ -164,50 +165,50 @@ def skew_LTL(A, overwrite_a=False, calc_L=True, calc_P=True):
     if calc_P:
         Pv = np.arange(n)
 
-    for k in range(n-2):
+    for k in range(n - 2):
         # First, find the largest entry in A[k+1:,k] and
         # permute it to A[k+1,k]
-        kp = k+1+np.abs(A[k+1:, k]).argmax()
+        kp = k + 1 + np.abs(A[k + 1 :, k]).argmax()
 
         # Check if we need to pivot
-        if kp != k+1:
+        if kp != k + 1:
             # interchange rows k+1 and kp
-            temp = A[k+1, k:].copy()
-            A[k+1, k:] = A[kp, k:]
+            temp = A[k + 1, k:].copy()
+            A[k + 1, k:] = A[kp, k:]
             A[kp, k:] = temp
 
             # Then interchange columns k+1 and kp
-            temp = A[k:, k+1].copy()
-            A[k:, k+1] = A[k:, kp]
+            temp = A[k:, k + 1].copy()
+            A[k:, k + 1] = A[k:, kp]
             A[k:, kp] = temp
 
             if calc_L:
                 # permute L accordingly
-                temp = L[k+1, 1:k+1].copy()
-                L[k+1, 1:k+1] = L[kp, 1:k+1]
-                L[kp, 1:k+1] = temp
+                temp = L[k + 1, 1 : k + 1].copy()
+                L[k + 1, 1 : k + 1] = L[kp, 1 : k + 1]
+                L[kp, 1 : k + 1] = temp
 
             if calc_P:
                 # accumulate the permutation matrix
-                temp = Pv[k+1]
-                Pv[k+1] = Pv[kp]
+                temp = Pv[k + 1]
+                Pv[k + 1] = Pv[kp]
                 Pv[kp] = temp
 
         # Now form the Gauss vector
-        if A[k+1, k] != 0.0:
-            tau = A[k+2:, k].copy()
-            tau /= A[k+1, k]
+        if A[k + 1, k] != 0.0:
+            tau = A[k + 2 :, k].copy()
+            tau /= A[k + 1, k]
 
             # clear eliminated row and column
-            A[k+2:, k] = 0.0
-            A[k, k+2:] = 0.0
+            A[k + 2 :, k] = 0.0
+            A[k, k + 2 :] = 0.0
 
             # Update the matrix block A(k+2:,k+2)
-            A[k+2:, k+2:] += np.outer(tau, A[k+2:, k+1])
-            A[k+2:, k+2:] -= np.outer(A[k+2:, k+1], tau)
+            A[k + 2 :, k + 2 :] += np.outer(tau, A[k + 2 :, k + 1])
+            A[k + 2 :, k + 2 :] -= np.outer(A[k + 2 :, k + 1], tau)
 
             if calc_L:
-                L[k+2:, k+1] = tau
+                L[k + 2 :, k + 1] = tau
 
     if calc_P:
         # form the permutation matrix as a sparse matrix
@@ -225,7 +226,7 @@ def skew_LTL(A, overwrite_a=False, calc_L=True, calc_P=True):
             return np.asmatrix(A)
 
 
-def pfaffian(A, overwrite_a=False, method='P', sign_only=False):
+def pfaffian(A, overwrite_a=False, method="P", sign_only=False):
     """ pfaffian(A, overwrite_a=False, method='P')
 
     Compute the Pfaffian of a real or complex skew-symmetric
@@ -237,12 +238,12 @@ def pfaffian(A, overwrite_a=False, method='P', sign_only=False):
     # Check if matrix is square
     assert A.shape[0] == A.shape[1] > 0
     # Check if it's skew-symmetric
-    assert abs((A+A.T).max()) < 1e-14, abs((A+A.T).max())
+    assert abs((A + A.T).max()) < 1e-14, abs((A + A.T).max())
     # Check that the method variable is appropriately set
-    assert method == 'P' or method == 'H'
-    if method == 'H' and sign_only:
+    assert method == "P" or method == "H"
+    if method == "H" and sign_only:
         raise Exception("Use `method='P'` when using `sign_only=True`")
-    if method == 'P':
+    if method == "P":
         return pfaffian_LTL(A, overwrite_a, sign_only)
     else:
         return pfaffian_householder(A, overwrite_a)
@@ -259,7 +260,7 @@ def pfaffian_LTL(A, overwrite_a=False, sign_only=False):
     # Check if matrix is square
     assert A.shape[0] == A.shape[1] > 0
     # Check if it's skew-symmetric
-    assert abs((A+A.T).max()) < 1e-14
+    assert abs((A + A.T).max()) < 1e-14
 
     n = A.shape[0]
     A = np.asarray(A)  # the slice views work only properly for arrays
@@ -273,40 +274,40 @@ def pfaffian_LTL(A, overwrite_a=False, sign_only=False):
 
     pfaffian_val = 1.0
 
-    for k in range(0, n-1, 2):
+    for k in range(0, n - 1, 2):
         # First, find the largest entry in A[k+1:,k] and
         # permute it to A[k+1,k]
-        kp = k+1+np.abs(A[k+1:, k]).argmax()
+        kp = k + 1 + np.abs(A[k + 1 :, k]).argmax()
 
         # Check if we need to pivot
-        if kp != k+1:
+        if kp != k + 1:
             # interchange rows k+1 and kp
-            temp = A[k+1, k:].copy()
-            A[k+1, k:] = A[kp, k:]
+            temp = A[k + 1, k:].copy()
+            A[k + 1, k:] = A[kp, k:]
             A[kp, k:] = temp
 
             # Then interchange columns k+1 and kp
-            temp = A[k:, k+1].copy()
-            A[k:, k+1] = A[k:, kp]
+            temp = A[k:, k + 1].copy()
+            A[k:, k + 1] = A[k:, kp]
             A[k:, kp] = temp
 
             # every interchange corresponds to a "-" in det(P)
             pfaffian_val *= -1
 
         # Now form the Gauss vector
-        if A[k+1, k] != 0.0:
-            tau = A[k, k+2:].copy()
-            tau /= A[k, k+1]
+        if A[k + 1, k] != 0.0:
+            tau = A[k, k + 2 :].copy()
+            tau /= A[k, k + 1]
 
             if sign_only:
-                pfaffian_val *= np.sign(A[k, k+1])
+                pfaffian_val *= np.sign(A[k, k + 1])
             else:
-                pfaffian_val *= A[k, k+1]
+                pfaffian_val *= A[k, k + 1]
 
-            if k+2 < n:
+            if k + 2 < n:
                 # Update the matrix block A(k+2:,k+2)
-                A[k+2:, k+2:] += np.outer(tau, A[k+2:, k+1])
-                A[k+2:, k+2:] -= np.outer(A[k+2:, k+1], tau)
+                A[k + 2 :, k + 2 :] += np.outer(tau, A[k + 2 :, k + 1])
+                A[k + 2 :, k + 2 :] -= np.outer(A[k + 2 :, k + 1], tau)
         else:
             # if we encounter a zero on the super/subdiagonal, the
             # Pfaffian is 0
@@ -331,7 +332,7 @@ def pfaffian_householder(A, overwrite_a=False):
     # Check if matrix is square
     assert A.shape[0] == A.shape[1] > 0
     # Check if it's skew-symmetric
-    assert abs((A+A.T).max()) < 1e-14
+    assert abs((A + A.T).max()) < 1e-14
 
     n = A.shape[0]
 
@@ -352,26 +353,26 @@ def pfaffian_householder(A, overwrite_a=False):
     if not overwrite_a:
         A = A.copy()
 
-    pfaffian_val = 1.
+    pfaffian_val = 1.0
 
-    for i in range(A.shape[0]-2):
+    for i in range(A.shape[0] - 2):
         # Find a Householder vector to eliminate the i-th column
-        v, tau, alpha = householder(A[i+1:, i])
-        A[i+1, i] = alpha
-        A[i, i+1] = -alpha
-        A[i+2:, i] = 0
-        A[i, i+2:] = 0
+        v, tau, alpha = householder(A[i + 1 :, i])
+        A[i + 1, i] = alpha
+        A[i, i + 1] = -alpha
+        A[i + 2 :, i] = 0
+        A[i, i + 2 :] = 0
 
         # Update the matrix block A(i+1:N,i+1:N)
-        w = tau * A[i+1:, i+1:] @ v.conj()
-        A[i+1:, i+1:] += np.outer(v, w) - np.outer(w, v)
+        w = tau * A[i + 1 :, i + 1 :] @ v.conj()
+        A[i + 1 :, i + 1 :] += np.outer(v, w) - np.outer(w, v)
 
         if tau != 0:
-            pfaffian_val *= 1-tau
+            pfaffian_val *= 1 - tau
         if i % 2 == 0:
             pfaffian_val *= -alpha
 
-    pfaffian_val *= A[n-2, n-1]
+    pfaffian_val *= A[n - 2, n - 1]
 
     return pfaffian_val
 
@@ -388,7 +389,8 @@ def pfaffian_schur(A, overwrite_a=False):
     """
 
     assert np.issubdtype(A.dtype, np.number) and not np.issubdtype(
-        A.dtype, np.complexfloating)
+        A.dtype, np.complexfloating
+    )
 
     assert A.shape[0] == A.shape[1] > 0
 
@@ -398,7 +400,7 @@ def pfaffian_schur(A, overwrite_a=False):
     if A.shape[0] % 2 == 1:
         return 0
 
-    (t, z) = la.schur(A, output='real', overwrite_a=overwrite_a)
+    (t, z) = la.schur(A, output="real", overwrite_a=overwrite_a)
     l = np.diag(t, 1)
     return np.prod(l[::2]) * la.det(z)
 
@@ -415,7 +417,7 @@ def pfaffian_sign(A, overwrite_a=False):
     # Check if matrix is square
     assert A.shape[0] == A.shape[1] > 0
     # Check if it's skew-symmetric
-    assert abs((A+A.T).max()) < 1e-14, abs((A+A.T).max())
+    assert abs((A + A.T).max()) < 1e-14, abs((A + A.T).max())
 
     return pfaffian_LTL_sign(A, overwrite_a)
 
@@ -431,7 +433,7 @@ def pfaffian_LTL_sign(A, overwrite_a=False):
     # Check if matrix is square
     assert A.shape[0] == A.shape[1] > 0
     # Check if it's skew-symmetric
-    assert abs((A+A.T).max()) < 1e-14
+    assert abs((A + A.T).max()) < 1e-14
 
     n = A.shape[0]
     A = np.asarray(A)  # the slice views work only properly for arrays
@@ -445,37 +447,37 @@ def pfaffian_LTL_sign(A, overwrite_a=False):
 
     pfaffian_val = 1.0
 
-    for k in range(0, n-1, 2):
+    for k in range(0, n - 1, 2):
         # First, find the largest entry in A[k+1:,k] and
         # permute it to A[k+1,k]
-        kp = k+1+np.abs(A[k+1:, k]).argmax()
+        kp = k + 1 + np.abs(A[k + 1 :, k]).argmax()
 
         # Check if we need to pivot
-        if kp != k+1:
+        if kp != k + 1:
             # interchange rows k+1 and kp
-            temp = A[k+1, k:].copy()
-            A[k+1, k:] = A[kp, k:]
+            temp = A[k + 1, k:].copy()
+            A[k + 1, k:] = A[kp, k:]
             A[kp, k:] = temp
 
             # Then interchange columns k+1 and kp
-            temp = A[k:, k+1].copy()
-            A[k:, k+1] = A[k:, kp]
+            temp = A[k:, k + 1].copy()
+            A[k:, k + 1] = A[k:, kp]
             A[k:, kp] = temp
 
             # every interchange corresponds to a "-" in det(P)
             pfaffian_val *= -1
 
         # Now form the Gauss vector
-        if A[k+1, k] != 0.0:
-            tau = A[k, k+2:].copy()
-            tau /= A[k, k+1]
+        if A[k + 1, k] != 0.0:
+            tau = A[k, k + 2 :].copy()
+            tau /= A[k, k + 1]
 
-            pfaffian_val *= A[k, k+1]
+            pfaffian_val *= A[k, k + 1]
 
-            if k+2 < n:
+            if k + 2 < n:
                 # Update the matrix block A(k+2:,k+2)
-                A[k+2:, k+2:] += np.outer(tau, A[k+2:, k+1])
-                A[k+2:, k+2:] -= np.outer(A[k+2:, k+1], tau)
+                A[k + 2 :, k + 2 :] += np.outer(tau, A[k + 2 :, k + 1])
+                A[k + 2 :, k + 2 :] -= np.outer(A[k + 2 :, k + 1], tau)
         else:
             # if we encounter a zero on the super/subdiagonal, the
             # Pfaffian is 0