diff --git a/gprMax/cmds_geometry/add_grass.py b/gprMax/cmds_geometry/add_grass.py
index 8689c9b2..04102295 100644
--- a/gprMax/cmds_geometry/add_grass.py
+++ b/gprMax/cmds_geometry/add_grass.py
@@ -111,7 +111,7 @@ class AddGrass(UserObjectGeometry):
             if ys == yf or zs == zf:
                 logger.exception(f"{self.__str__()} dimensions are not specified correctly")
                 raise ValueError
-            if xs != volume.xs and xs != volume.xf:
+            if xs not in [volume.xs and volume.xf]:
                 logger.exception(f"{self.__str__()} must specify external surfaces on a fractal box")
                 raise ValueError
             fractalrange = (round_value(limits[0] / grid.dx), round_value(limits[1] / grid.dx))
@@ -133,10 +133,10 @@ class AddGrass(UserObjectGeometry):
                 requestedsurface = "xplus"
 
         elif ys == yf:
-            if xs == xf or zs == zf:
+            if zs == zf:
                 logger.exception(f"{self.__str__()} dimensions are not specified correctly")
                 raise ValueError
-            if ys != volume.ys and ys != volume.yf:
+            if ys not in [volume.ys and volume.yf]:
                 logger.exception(f"{self.__str__()} must specify external surfaces on a fractal box")
                 raise ValueError
             fractalrange = (round_value(limits[0] / grid.dy), round_value(limits[1] / grid.dy))
@@ -158,10 +158,7 @@ class AddGrass(UserObjectGeometry):
                 requestedsurface = "yplus"
 
         elif zs == zf:
-            if xs == xf or ys == yf:
-                logger.exception(f"{self.__str__()} dimensions are not specified correctly")
-                raise ValueError
-            if zs != volume.zs and zs != volume.zf:
+            if zs not in [volume.zs and volume.zf]:
                 logger.exception(f"{self.__str__()} must specify external surfaces on a fractal box")
                 raise ValueError
             fractalrange = (round_value(limits[0] / grid.dz), round_value(limits[1] / grid.dz))
@@ -238,7 +235,7 @@ class AddGrass(UserObjectGeometry):
         surface.grass.append(g)
 
         # Check to see if grass has been already defined as a material
-        if not any(x.ID == "grass" for x in grid.materials):
+        if all(x.ID == "grass" for x in grid.materials):
             create_grass(grid)
 
         # Check if time step for model is suitable for using grass
diff --git a/gprMax/cmds_geometry/add_surface_roughness.py b/gprMax/cmds_geometry/add_surface_roughness.py
index bb6f9e36..00eb86cf 100644
--- a/gprMax/cmds_geometry/add_surface_roughness.py
+++ b/gprMax/cmds_geometry/add_surface_roughness.py
@@ -123,7 +123,7 @@ class AddSurfaceRoughness(UserObjectGeometry):
             if ys == yf or zs == zf:
                 logger.exception(f"{self.__str__()} dimensions are not specified correctly")
                 raise ValueError
-            if xs != volume.xs and xs != volume.xf:
+            if xs not in [volume.xs, volume.xf]:
                 logger.exception(f"{self.__str__()} can only be used on the external " + "surfaces of a fractal box")
                 raise ValueError
             fractalrange = (round_value(limits[0] / grid.dx), round_value(limits[1] / grid.dx))
@@ -151,10 +151,10 @@ class AddSurfaceRoughness(UserObjectGeometry):
                 requestedsurface = "xplus"
 
         elif ys == yf:
-            if xs == xf or zs == zf:
+            if zs == zf:
                 logger.exception(f"{self.__str__()} dimensions are not specified correctly")
                 raise ValueError
-            if ys != volume.ys and ys != volume.yf:
+            if ys not in [volume.ys and volume.yf]:
                 logger.exception(f"{self.__str__()} can only be used on the external " + "surfaces of a fractal box")
                 raise ValueError
             fractalrange = (round_value(limits[0] / grid.dy), round_value(limits[1] / grid.dy))
@@ -182,10 +182,7 @@ class AddSurfaceRoughness(UserObjectGeometry):
                 requestedsurface = "yplus"
 
         elif zs == zf:
-            if xs == xf or ys == yf:
-                logger.exception(f"{self.__str__()} dimensions are not specified correctly")
-                raise ValueError
-            if zs != volume.zs and zs != volume.zf:
+            if zs not in [volume.zs and volume.zf]:
                 logger.exception(f"{self.__str__()} can only be used on the external " + "surfaces of a fractal box")
                 raise ValueError
             fractalrange = (round_value(limits[0] / grid.dz), round_value(limits[1] / grid.dz))
diff --git a/gprMax/cmds_geometry/fractal_box.py b/gprMax/cmds_geometry/fractal_box.py
index a205415f..1aea4a8c 100644
--- a/gprMax/cmds_geometry/fractal_box.py
+++ b/gprMax/cmds_geometry/fractal_box.py
@@ -144,14 +144,13 @@ class FractalBox(UserObjectGeometry):
                     f"{self.__str__()} must be used with more than " + "one material from the mixing model."
                 )
                 raise ValueError
-            if isinstance(mixingmodel, ListMaterial):
-                if nbins > len(mixingmodel.mat):
-                    logger.exception(
-                        f"{self.__str__()} too many materials/bins "
-                        + "requested compared to materials in "
-                        + "mixing model."
-                    )
-                    raise ValueError
+            if isinstance(mixingmodel, ListMaterial) and nbins > len(mixingmodel.mat):
+                logger.exception(
+                    f"{self.__str__()} too many materials/bins "
+                    + "requested compared to materials in "
+                    + "mixing model."
+                )
+                raise ValueError
             # Create materials from mixing model as number of bins now known
             # from fractal_box command.
             mixingmodel.calculate_properties(nbins, grid)
diff --git a/gprMax/cmds_multiuse.py b/gprMax/cmds_multiuse.py
index b5127cb6..0ba36075 100644
--- a/gprMax/cmds_multiuse.py
+++ b/gprMax/cmds_multiuse.py
@@ -385,7 +385,7 @@ class HertzianDipole(UserObjectMulti):
         p2 = uip.round_to_grid_static_point(p1)
 
         # Check if there is a waveformID in the waveforms list
-        if not any(x.ID == waveform_id for x in grid.waveforms):
+        if all(x.ID == waveform_id for x in grid.waveforms):
             logger.exception(f"{self.params_str()} there is no waveform " + f"with the identifier {waveform_id}.")
             raise ValueError
 
@@ -526,7 +526,7 @@ class MagneticDipole(UserObjectMulti):
         p2 = uip.round_to_grid_static_point(p1)
 
         # Check if there is a waveformID in the waveforms list
-        if not any(x.ID == waveform_id for x in grid.waveforms):
+        if all(x.ID == waveform_id for x in grid.waveforms):
             logger.exception(f"{self.params_str()} there is no waveform " + f"with the identifier {waveform_id}.")
             raise ValueError
 
diff --git a/gprMax/config.py b/gprMax/config.py
index d75844c1..fbbb7628 100644
--- a/gprMax/config.py
+++ b/gprMax/config.py
@@ -73,7 +73,7 @@ class ModelConfig:
         #     N.B. This will happen if the requested snapshots are too large to
         #           fit on the memory of the GPU. If True this will slow
         #           performance significantly.
-        if sim_config.general["solver"] == "cuda" or sim_config.general["solver"] == "opencl":
+        if sim_config.general["solver"] in ["cuda", "opencl"]:
             if sim_config.general["solver"] == "cuda":
                 devs = sim_config.args.gpu
             elif sim_config.general["solver"] == "opencl":
diff --git a/gprMax/contexts.py b/gprMax/contexts.py
index 38a10856..4f3e66a0 100644
--- a/gprMax/contexts.py
+++ b/gprMax/contexts.py
@@ -51,7 +51,6 @@ class Context:
             results: dict that can contain useful results/data from simulation.
         """
 
-        results = {}
         self.tsimstart = timer()
         self.print_logo_copyright()
         print_host_info(config.sim_config.hostinfo)
@@ -86,7 +85,7 @@ class Context:
         self.tsimend = timer()
         self.print_sim_time_taken()
 
-        return results
+        return {}
 
     def print_logo_copyright(self):
         """Prints gprMax logo, version, and copyright/licencing information."""
diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py
index 9b96b14c..ba61add0 100644
--- a/gprMax/gprMax.py
+++ b/gprMax/gprMax.py
@@ -217,10 +217,9 @@ def run_main(args):
     # MPI running with (OpenMP/CUDA/OpenCL)
     if config.sim_config.args.mpi:
         context = MPIContext()
-        results = context.run()
     # Standard running (OpenMP/CUDA/OpenCL)
     else:
         context = Context()
-        results = context.run()
 
+    results = context.run()
     return results
diff --git a/gprMax/grid.py b/gprMax/grid.py
index 4da21882..e0744e87 100644
--- a/gprMax/grid.py
+++ b/gprMax/grid.py
@@ -417,7 +417,7 @@ def dispersion_analysis(G):
     # Find maximum significant frequency
     if G.waveforms:
         for waveform in G.waveforms:
-            if waveform.type == "sine" or waveform.type == "contsine":
+            if waveform.type in ["sine", "contsine"]:
                 results["maxfreq"].append(4 * waveform.freq)
 
             elif waveform.type == "impulse":
diff --git a/gprMax/hash_cmds_file.py b/gprMax/hash_cmds_file.py
index 66e26dad..d8529763 100644
--- a/gprMax/hash_cmds_file.py
+++ b/gprMax/hash_cmds_file.py
@@ -322,14 +322,13 @@ def check_cmd_names(processedlines, checkessential=True):
 
         lindex += 1
 
-    if checkessential:
-        if countessentialcmds < len(essentialcmds):
-            logger.exception(
-                "Your input file is missing essential commands "
-                + "required to run a model. Essential commands are: "
-                + ", ".join(essentialcmds)
-            )
-            raise SyntaxError
+    if checkessential and countessentialcmds < len(essentialcmds):
+        logger.exception(
+            "Your input file is missing essential commands "
+            + "required to run a model. Essential commands are: "
+            + ", ".join(essentialcmds)
+        )
+        raise SyntaxError
 
     return singlecmds, multiplecmds, geometry
 
diff --git a/gprMax/updates.py b/gprMax/updates.py
index 34973420..9adece69 100644
--- a/gprMax/updates.py
+++ b/gprMax/updates.py
@@ -1209,7 +1209,6 @@ class OpenCLUpdates:
                 self.grid.Hy_dev,
                 self.grid.Hz_dev,
             )
-            event.wait()
 
         # If there are any dispersive materials do 1st part of dispersive update
         # (it is split into two parts as it requires present and updated electric field values).
@@ -1231,7 +1230,8 @@ class OpenCLUpdates:
                 self.grid.Hy_dev,
                 self.grid.Hz_dev,
             )
-            event.wait()
+
+        event.wait()
 
     def update_electric_pml(self):
         """Updates electric field components with the PML correction."""
@@ -1319,7 +1319,6 @@ class OpenCLUpdates:
         # if iteration == self.grid.iterations - 1:
         #     return self.drv.mem_get_info()[1] - self.drv.mem_get_info()[0]
         logger.debug("Look at memory estimate for pyopencl")
-        pass
 
     def calculate_solve_time(self):
         """Calculates solving time for model."""
diff --git a/testing/benchmarking/bench_simple.py b/testing/benchmarking/bench_simple.py
index 7fec7406..1e7d1edb 100644
--- a/testing/benchmarking/bench_simple.py
+++ b/testing/benchmarking/bench_simple.py
@@ -3,9 +3,12 @@
     receiver at the centre.
 """
 
+
 from pathlib import Path
 import gprMax
 
+import itertools
+
 # File path for output
 fn = Path(__file__)
 
@@ -17,35 +20,34 @@ ompthreads = [1, 2, 4, 8, 16, 32, 64, 128]
 
 scenes = []
 
-for d in domains:
-    for threads in ompthreads:
-        # Discretisation
-        dl = 0.001
+# Discretisation
+dl = 0.001
 
-        # Domain
-        x = d
-        y = x
-        z = x
+for d, threads in itertools.product(domains, ompthreads):
+    # Domain
+    x = d
+    y = x
+    z = x
 
-        scene = gprMax.Scene()
+    scene = gprMax.Scene()
 
-        title = gprMax.Title(name=fn.with_suffix("").name)
-        domain = gprMax.Domain(p1=(x, y, z))
-        dxdydz = gprMax.Discretisation(p1=(dl, dl, dl))
-        time_window = gprMax.TimeWindow(time=3e-9)
-        wv = gprMax.Waveform(wave_type="gaussiandotnorm", amp=1, freq=900e6, id="MySource")
-        src = gprMax.HertzianDipole(p1=(x / 2, y / 2, z / 2), polarisation="x", waveform_id="MySource")
+    title = gprMax.Title(name=fn.with_suffix("").name)
+    domain = gprMax.Domain(p1=(x, y, z))
+    dxdydz = gprMax.Discretisation(p1=(dl, dl, dl))
+    time_window = gprMax.TimeWindow(time=3e-9)
+    wv = gprMax.Waveform(wave_type="gaussiandotnorm", amp=1, freq=900e6, id="MySource")
+    src = gprMax.HertzianDipole(p1=(x / 2, y / 2, z / 2), polarisation="x", waveform_id="MySource")
 
-        omp = gprMax.OMPThreads(n=threads)
+    omp = gprMax.OMPThreads(n=threads)
 
-        scene.add(title)
-        scene.add(domain)
-        scene.add(dxdydz)
-        scene.add(time_window)
-        scene.add(wv)
-        scene.add(src)
-        scene.add(omp)
-        scenes.append(scene)
+    scene.add(title)
+    scene.add(domain)
+    scene.add(dxdydz)
+    scene.add(time_window)
+    scene.add(wv)
+    scene.add(src)
+    scene.add(omp)
+    scenes.append(scene)
 
 # Run model
 gprMax.run(scenes=scenes, n=len(scenes), geometry_only=False, outputfile=fn, gpu=None)
diff --git a/toolboxes/AntennaPatterns/plot_fields.py b/toolboxes/AntennaPatterns/plot_fields.py
index 1d14b6de..b943369a 100644
--- a/toolboxes/AntennaPatterns/plot_fields.py
+++ b/toolboxes/AntennaPatterns/plot_fields.py
@@ -61,9 +61,9 @@ if epsr:
     wavelength = v1 / f
 
 # Print some useful information
-logger.info("Centre frequency: {} GHz".format(f / 1e9))
+logger.info(f"Centre frequency: {f / 1000000000.0} GHz")
 if epsr:
-    logger.info("Critical angle for Er {} is {} degrees".format(epsr, thetac))
+    logger.info(f"Critical angle for Er {epsr} is {thetac} degrees")
     logger.info("Wavelength: {:.3f} m".format(wavelength))
     logger.info(
         "Observation distance(s) from {:.3f} m ({:.1f} wavelengths) to {:.3f} m ({:.1f} wavelengths)".format(
@@ -139,7 +139,7 @@ leg = ax.legend(
 [legobj.set_linewidth(2) for legobj in leg.legendHandles]
 
 # Save a pdf of the plot
-savename = os.path.splitext(args.numpyfile)[0] + ".pdf"
+savename = f"{os.path.splitext(args.numpyfile)[0]}.pdf"
 fig.savefig(savename, dpi=None, format="pdf", bbox_inches="tight", pad_inches=0.1)
 # savename = os.path.splitext(args.numpyfile)[0] + '.png'
 # fig.savefig(savename, dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
diff --git a/toolboxes/DebyeFit/Debye_Fit.py b/toolboxes/DebyeFit/Debye_Fit.py
index 28faf219..221ca1ac 100644
--- a/toolboxes/DebyeFit/Debye_Fit.py
+++ b/toolboxes/DebyeFit/Debye_Fit.py
@@ -123,9 +123,7 @@ class Relaxation(object):
         """
         print(f"Approximating {self.name}" f" using {self.number_of_debye_poles} Debye poles")
         print(f"{self.name} parameters: ")
-        s = ""
-        for k, v in self.params.items():
-            s += f"{k:10s} = {v}\n"
+        s = "".join(f"{k:10s} = {v}\n" for k, v in self.params.items())
         print(s)
         return f"{self.name}:\n{s}"
 
@@ -232,10 +230,10 @@ class Relaxation(object):
             "#material: {} {} {} {} {}\n".format(ee, self.sigma, self.mu, self.mu_sigma, self.material_name)
         )
         print(material_prop[0], end="")
-        dispersion_prop = "#add_dispersion_debye: {}".format(len(tau))
+        dispersion_prop = f"#add_dispersion_debye: {len(tau)}"
         for i in range(len(tau)):
-            dispersion_prop += " {} {}".format(weights[i], 10 ** tau[i])
-        dispersion_prop += " {}".format(self.material_name)
+            dispersion_prop += f" {weights[i]} {10**tau[i]}"
+        dispersion_prop += f" {self.material_name}"
         print(dispersion_prop)
         material_prop.append(dispersion_prop + "\n")
         return material_prop
@@ -312,11 +310,10 @@ class Relaxation(object):
             file_path = os.path.join("user_libs", "materials", "my_materials.txt")
         else:
             sys.exit("Cannot save material properties " f"in {os.path.join(fdir, 'my_materials.txt')}!")
-        fileH = open(file_path, "a")
-        fileH.write(f"## {output[0].split(' ')[-1]}")
-        fileH.writelines(output)
-        fileH.write("\n")
-        fileH.close()
+        with open(file_path, "a") as fileH:
+            fileH.write(f"## {output[0].split(' ')[-1]}")
+            fileH.writelines(output)
+            fileH.write("\n")
         print(f"Material properties save at: {file_path}")
 
 
@@ -613,7 +610,7 @@ class Crim(Relaxation):
         print(f"Approximating Complex Refractive Index Model (CRIM)" f" using {self.number_of_debye_poles} Debye poles")
         print("CRIM parameters: ")
         for i in range(len(self.volumetric_fractions)):
-            print("Material {}.:".format(i + 1))
+            print(f"Material {i + 1}.:")
             print("---------------------------------")
             print(f"{'Vol. fraction':>27s} = {self.volumetric_fractions[i]}")
             print(f"{'e_inf':>27s} = {self.materials[i][0]}")
diff --git a/toolboxes/DebyeFit/optimization.py b/toolboxes/DebyeFit/optimization.py
index 342d880d..ff94b2f6 100644
--- a/toolboxes/DebyeFit/optimization.py
+++ b/toolboxes/DebyeFit/optimization.py
@@ -474,7 +474,6 @@ def DLS(logt, rl, im, freq):
     rp, ip = np.matmul(d.real, x[np.newaxis].T).T[0], np.matmul(d.imag, x[np.newaxis].T).T[0]
     cost_i = np.sum(np.abs(ip - im)) / len(im)
     ee = np.mean(rl - rp)
-    if ee < 1:
-        ee = 1
+    ee = max(ee, 1)
     cost_r = np.sum(np.abs(rp + ee - rl)) / len(im)
     return cost_i, cost_r, x, ee, rp, ip
diff --git a/toolboxes/Plotting/plot_Ascan.py b/toolboxes/Plotting/plot_Ascan.py
index 3990e0f0..9af059cd 100644
--- a/toolboxes/Plotting/plot_Ascan.py
+++ b/toolboxes/Plotting/plot_Ascan.py
@@ -79,10 +79,9 @@ def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False, save=False):
         time = np.linspace(0, (iterations - 1) * dt, num=iterations)
 
         # Check for single output component when doing a FFT
-        if fft:
-            if not len(outputs) == 1:
-                logger.exception("A single output must be specified when using " + "the -fft option")
-                raise ValueError
+        if fft and not len(outputs) == 1:
+            logger.exception("A single output must be specified when using " + "the -fft option")
+            raise ValueError
 
         # New plot for each receiver
         for rx in range(1, nrx + 1):
diff --git a/toolboxes/Plotting/plot_source_wave.py b/toolboxes/Plotting/plot_source_wave.py
index cd838143..fa86e649 100644
--- a/toolboxes/Plotting/plot_source_wave.py
+++ b/toolboxes/Plotting/plot_source_wave.py
@@ -86,19 +86,13 @@ def mpl_plot(w, timewindow, dt, iterations, fft=False, save=False):
     logging.info(f"Type: {w.type}")
     logging.info(f"Maximum (absolute) amplitude: {np.max(np.abs(waveform)):g}")
 
-    if w.freq and not w.type == "gaussian" and not w.type == "impulse":
+    if w.freq and w.type != "gaussian" and w.type != "impulse":
         logging.info(f"Centre frequency: {w.freq:g} Hz")
 
-    if (
-        w.type == "gaussian"
-        or w.type == "gaussiandot"
-        or w.type == "gaussiandotnorm"
-        or w.type == "gaussianprime"
-        or w.type == "gaussiandoubleprime"
-    ):
+    if w.type in ["gaussian", "gaussiandot", "gaussiandotnorm", "gaussianprime", "gaussiandoubleprime"]:
         delay = 1 / w.freq
         logging.info(f"Time to centre of pulse: {delay:g} s")
-    elif w.type == "gaussiandotdot" or w.type == "gaussiandotdotnorm" or w.type == "ricker":
+    elif w.type in ["gaussiandotdot", "gaussiandotdotnorm", "ricker"]:
         delay = np.sqrt(2) / w.freq
         logging.info(f"Time to centre of pulse: {delay:g} s")
 
diff --git a/toolboxes/STLtoVoxel/perimeter.py b/toolboxes/STLtoVoxel/perimeter.py
index ea28a047..5952c1ec 100644
--- a/toolboxes/STLtoVoxel/perimeter.py
+++ b/toolboxes/STLtoVoxel/perimeter.py
@@ -36,12 +36,9 @@ def generate_y(p1, p2, x):
 
 def paint_y_axis(lines, pixels, x):
     is_black = False
-    target_ys = list(map(lambda line: int(generate_y(line[0], line[1], x)), lines))
-    target_ys.sort()
+    target_ys = sorted(map(lambda line: int(generate_y(line[0], line[1], x)), lines))
     if len(target_ys) % 2:
-        distances = []
-        for i in range(len(target_ys) - 1):
-            distances.append(target_ys[i + 1] - target_ys[i])
+        distances = [target_ys[i + 1] - target_ys[i] for i in range(len(target_ys) - 1)]
         # https://stackoverflow.com/a/17952763
         min_idx = -min((x, -i) for i, x in enumerate(distances))[1]
         del target_ys[min_idx]
@@ -54,7 +51,7 @@ def paint_y_axis(lines, pixels, x):
         pixels[target_y][x] = True
         is_black = not is_black
         yi = target_y
-    assert is_black is False, "an error has occured at x%s" % x
+    assert is_black is False, f"an error has occured at x{x}"
 
 
 def generate_line_events(line_list):
diff --git a/toolboxes/STLtoVoxel/slice.py b/toolboxes/STLtoVoxel/slice.py
index 7767f826..32802c89 100644
--- a/toolboxes/STLtoVoxel/slice.py
+++ b/toolboxes/STLtoVoxel/slice.py
@@ -1,3 +1,4 @@
+import itertools
 import multiprocessing as mp
 import sys
 
@@ -90,10 +91,7 @@ def triangle_to_intersecting_lines(triangle, height, pixels, lines):
             y = int(same[0][1])
             pixels[y][x] = True
     else:
-        cross_lines = []
-        for a in above:
-            for b in below:
-                cross_lines.append((b, a))
+        cross_lines = [(b, a) for a, b in itertools.product(above, below)]
         side1 = where_line_crosses_z(cross_lines[0][0], cross_lines[0][1], height)
         side2 = where_line_crosses_z(cross_lines[1][0], cross_lines[1][1], height)
         lines.append((side1, side2))