diff --git a/.gitignore b/.gitignore
index 5240c39e..41c58297 100644
--- a/.gitignore
+++ b/.gitignore
@@ -50,4 +50,5 @@ DVASP_MACE_comparison/
 pca_model*.pkl
 .benchmarks/
 .local-pre-commit-config.yaml
-*.tgz
\ No newline at end of file
+*.tgz
+.ipynb_checkpoints
\ No newline at end of file
diff --git a/app/inputs.f90 b/app/inputs.f90
index 6bc6be7b..09dc9833 100644
--- a/app/inputs.f90
+++ b/app/inputs.f90
@@ -126,7 +126,7 @@ subroutine set_global_vars()
              write(0,'("&
                   &ERROR: No input filename supplied, &
                   &but the flag ''-f'' was used&
-                  &")')
+             &")')
              infilename_do: do j = 1, 3
                 write(6,'("Please supply an input filename:")')
                 read(5,'(A)') input_file
@@ -157,7 +157,7 @@ subroutine set_global_vars()
           write(6,'("-----------------FILE-NAME-FLAGS-----------------")')
           write(6,'(2X,"-f<STR>         : Input structure file name (&
                &Default = (empty)&
-               &).")')
+          &).")')
           stop
        end if
     end do flagloop
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 9b89486a..2776a729 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -15,8 +15,10 @@ An example
 
 .. code-block:: python
 
-    # A simple example of how to use RAFFLE to generate 10 structures of diamond
+    # A simple example of how to use RAFFLE to generate 10 structures of diamond and write them to a single file
     from ase import Atoms
+    from ase.io import write
+    from ase.calculators.singlepoint import SinglePointCalculator
     from raffle.generator import raffle_generator
     from mace.calculators import mace_mp
 
@@ -41,6 +43,16 @@ An example
         generator.distributions.update(structures[num_structures_old:])
         num_structures_old = len(structures)
 
+    structures = generator.get_structures(calc)
+    for structure in structures:
+        structure.calc = SinglePointCalculator(
+            structure,
+            energy=structure.get_potential_energy(),
+            forces=structure.get_forces()
+        )
+ 
+    write('structures.traj', structures)
+
 .. toctree::
    :maxdepth: 3
    :caption: Contents:
diff --git a/docs/source/tutorials/BaTiO3_tutorial.rst b/docs/source/tutorials/BaTiO3_tutorial.rst
deleted file mode 100644
index 3c136a68..00000000
--- a/docs/source/tutorials/BaTiO3_tutorial.rst
+++ /dev/null
@@ -1,5 +0,0 @@
-.. BaTiO3:
-
-========================
-BaTiO\ :sub:`3` tutorial
-========================
diff --git a/docs/source/tutorials/C-MgO_tutorial.rst b/docs/source/tutorials/C-MgO_tutorial.rst
deleted file mode 100644
index 7a9b11e3..00000000
--- a/docs/source/tutorials/C-MgO_tutorial.rst
+++ /dev/null
@@ -1,5 +0,0 @@
-.. C-MgO:
-
-=========================
-Graphene-encapsulated MgO
-=========================
\ No newline at end of file
diff --git a/docs/source/tutorials/Si-Ge_tutorial.rst b/docs/source/tutorials/Si-Ge_tutorial.rst
new file mode 100644
index 00000000..2da2b308
--- /dev/null
+++ b/docs/source/tutorials/Si-Ge_tutorial.rst
@@ -0,0 +1,116 @@
+.. si-ge:
+
+========================
+Si|Ge interface tutorial
+========================
+
+This tutorial will guide you through the process of performing RAFFLE-based structure search for an Si|Ge interface.
+
+The tutorial is designed to show how a RAFFLE generator given no prior knowledge can learn bonding and statistically available interface configurations.
+Statistically available configurations are those that are likely to exist due to energetic values close to the ground state, within thermal and entropy conditions.
+
+The example script can be found in the following directory:
+
+.. code-block:: bash
+
+    raffle/example/python_pkg/Si-Ge_learn/learn.py
+
+We recommend reading through the file and running it to understand the process of learning and generating structures.
+However, we will provide a brief overview of the script here.
+
+First, we must import the required packages:
+
+.. code-block:: python
+
+    from ase import Atoms
+    from raffle.generator import raffle_generator
+    from mace.calculators import mace_mp
+    import numpy as np
+
+Next, we need to set up the RAFFLE generator and the calculator to calculate the energies of the structures.
+In this example, we use the CHGNet calculator:
+
+.. code-block:: python
+
+    generator = raffle_generator()
+
+    calc = mace_mp(model="mace-mpa-0-medium.model")
+
+Note, choice of calculator is important.
+The calculator should be valid within and near the chemical environment of the structures being generated; in this case, Si and Ge bonding.
+If this is the case, it can accurately identify local minima and provide a good representation of the energy landscape.
+For this example, we use the [MACE-MPA-0 calculator](https://github.com/ACEsuit/mace-mp/releases/tag/mace_mpa_0), which, from preliminary testing, has been found to be suitable for Si and Ge bonding.
+To use the MACE-MPA-0 calculator, you will need to download the model file from the link provided and place it in the same directory as the script (and install the MACE package using `pip install mace-torch`, version 0.3.10 or later).
+The CHGNet calculator was not found to be suitable for Si and Ge interface bonding.
+
+Then, we need to create the host structure.
+Here, an abrupt Si|Ge interface is generated.
+This is the base structure that will be added to in order to generate the structures.
+A vacuum region of 5.54 Å is set up between the Si and Ge regions, in which the atoms will be placed by RAFFLE.
+
+.. code-block:: python
+
+    Si_bulk = build.bulk("Si", crystalstructure="diamond", a=5.43)
+    Si_cubic = build.make_supercell(Si_bulk, [[-1, 1, 1], [1, -1, 1], [1, 1, -1]])
+    Ge_bulk = build.bulk("Ge", crystalstructure="diamond", a=5.65)
+    Ge_cubic = build.make_supercell(Ge_bulk, [[-1, 1, 1], [1, -1, 1], [1, 1, -1]])
+
+    Si_supercell = build.make_supercell(Si_cubic, [[2, 0, 0], [0, 2, 0], [0, 0, 1]])
+    Ge_supercell = build.make_supercell(Ge_cubic, [[2, 0, 0], [0, 2, 0], [0, 0, 1]])
+
+    Si_surface = build.surface(Si_supercell, indices=(0, 0, 1), layers=2)
+    Ge_surface = build.surface(Ge_supercell, indices=(0, 0, 1), layers=2)
+
+    host = build.stack(Si_surface, Ge_surface, axis=2, distance= 5.43/2 + 5.65/2)
+    cell[2, 2] -= 3.8865
+    host.set_cell(cell, scale_atoms=False)
+
+
+The script then sets parameters for the generator and provides an initial database.
+Note, this database only contains prior information of Si-Si and Ge-Ge bonding, and no information about Si-Ge bonding.
+This is to demonstrate the ability of RAFFLE to learn from scratch.
+
+.. code-block:: python
+
+    Si_bulk.calc = calc
+    Ge_bulk.calc = calc
+    generator.distributions.set_element_energies(
+        {
+            'Si': Si_bulk.get_potential_energy() / len(Si_bulk),
+            'Ge': Ge_bulk.get_potential_energy() / len(Ge_bulk),
+        }
+    )
+
+    # set energy scale
+    generator.distributions.set_kBT(0.2)
+
+    # set the distribution function widths (2-body, 3-body, 4-body)
+    generator.distributions.set_width([0.04, np.pi/160.0, np.pi/160.0])
+
+    # set the initial database
+    initial_database = [Si_bulk, Ge_bulk]
+    generator.distributions.create(initial_database)
+
+Finally, the script generates structures using the generator.
+The generator is given the host structures.
+Finally, the generator is run for each host structure, providing a unique stoichiometry each time and using a custom method ratio.
+The 
+
+.. code-block:: python
+
+    generator.set_host(host)
+    generator.set_bounds([[0, 0, 0.34], [1, 1, 0.52]])
+    for iter in range(40):
+        # generate the structures
+        structures, exit_code = generator.generate(
+            num_structures = 5,
+            stoichiometry = { 'Si': 16, 'Ge': 16 },
+            seed = iter,
+            method_ratio = {"void": 0.1, "rand": 0.01, "walk": 0.25, "grow": 0.25, "min": 1.0},
+            verbose = 0,
+            calc = calc
+        )
+        generator.distributions.update(structures)
+
+    structures = generator.get_structures()
+    write('structures.traj', structures)
diff --git a/docs/source/tutorials/aluminium_tutorial.rst b/docs/source/tutorials/aluminium_tutorial.rst
new file mode 100644
index 00000000..cea6c158
--- /dev/null
+++ b/docs/source/tutorials/aluminium_tutorial.rst
@@ -0,0 +1,86 @@
+.. aluminium:
+
+==================
+Aluminium tutorial
+==================
+
+This tutorial will guide you through the process of performing RAFFLE-based structure search for the bulk phases of aluminium.
+
+The tutorial is designed to show how a RAFFLE generator given no prior knowledge can learn bonding and identify known phases.
+This is not an expected use-case of RAFFLE due to its application to a bulk system, but still demonstrates the expected workflow and capabilities.
+
+The example script can be found in the following directory:
+
+.. code-block:: bash
+
+    raffle/example/python_pkg/Al_learn/learn.py
+
+We recommend reading through the file and running it to understand the process of learning and generating structures.
+However, we will provide a brief overview of the script here.
+
+First, we must import the required packages:
+
+.. code-block:: python
+
+    from ase import Atoms
+    from raffle.generator import raffle_generator
+    from chgnet.model import CHGNetCalculator
+    import numpy as np
+
+Next, we need to set up the RAFFLE generator and the calculator to calculate the energies of the structures.
+In this example, we use the CHGNet calculator:
+
+.. code-block:: python
+
+    generator = raffle_generator()
+
+    calc = CHGNetCalculator()
+
+Then, we need to create the host structure.
+Here, a set of host cells are generated to represent the multiple potential bulk phases of aluminium.
+These are then used to generate the structures.
+
+.. code-block:: python
+
+    crystal_structures = ['orthorhombic', 'hcp']
+    hosts = []
+    for crystal_structure in crystal_structures:
+        for a in np.linspace(3.1, 5.4, num=6):
+            atom = build.bulk(
+                    name = 'Al',
+                    crystalstructure = crystal_structure,
+                    a = a, b = a, c = a,
+            )
+            hosts.append(Atoms(
+                'Al',
+                positions = [(0, 0, 0)],
+                cell = atom.get_cell(),
+                pbc = True,
+                calculator = calc
+            ))
+
+The script then sets parameters for the generator and provides an initial database.
+Note, this database is effectively empty, as it only contains a single structure, which is an isolated aluminium atom.
+
+.. code-block:: python
+
+    initial_database = [Atoms('Al', positions=[(0, 0, 0)], cell=[8, 8, 8], pbc=True)]
+    initial_database[0].calc = calc
+    generator.distributions.create(initial_database)
+
+Finally, the script generates structures using the generator.
+The generator is given the host structures.
+Finally, the generator is run for each host structure, providing a unique stoichiometry each time and using a custom method ratio.
+
+.. code-block:: python
+
+    for host in hosts:
+        generator.set_host(host)
+        generator.generate(
+            num_structures = 5,
+            stoichiometry = { 'Al': num_atoms },
+            method_ratio = {"void": 0.5, "rand": 0.001, "walk": 0.5, "grow": 0.0, "min": 1.0},
+        )
+
+    structures = generator.get_structures()
+    write('structures.traj', structures)
diff --git a/docs/source/tutorials/diamond_tutorial.rst b/docs/source/tutorials/diamond_tutorial.rst
index ef5410c0..5420e90d 100644
--- a/docs/source/tutorials/diamond_tutorial.rst
+++ b/docs/source/tutorials/diamond_tutorial.rst
@@ -61,7 +61,6 @@ This is the base structure that will be added to in order to generate the struct
                 3.5607451090903233, 3.5607451090903233, 7.1214902182
             ], pbc=True
         )
-    )
 
     host.calc = calc
     generator.set_host(host)
@@ -121,12 +120,12 @@ We are now ready to generate structures using the database of structures.
 .. code-block:: python
 
     num_structures_old = 0
-    generator.generate(
+    structures, exit_code = generator.generate(
         num_structures = 1,
         stoichiometry = { 'C': 8 },
         method_ratio = {"void":0.0001, "min":1.0},
+        calc = calc
     )
-    structures = generator.get_structures(calc)
 
 We should now have a structure of diamond.
 This structure can be visualised using the ASE package.
diff --git a/docs/source/tutorials/graphite_tutorial.rst b/docs/source/tutorials/graphite_tutorial.rst
index 0a2e6aa4..eb95f92a 100644
--- a/docs/source/tutorials/graphite_tutorial.rst
+++ b/docs/source/tutorials/graphite_tutorial.rst
@@ -2,4 +2,128 @@
 
 =================
 Graphite tutorial
-=================
\ No newline at end of file
+=================
+
+This tutorial will guide you through the process of reconstructing graphite from a defected graphite cell (one missing layer).
+This tutorial follows the same structure as the `Diamond tutorial <diamond>`_.
+
+The tutorial is designed to show how RAFFLE learns from a database and uses this information to generate structures.
+This is not an expected use-case of RAFFLE, but merely a demonstration of its capabilities to rebuild a structure from a defected cell.
+
+The example files can be found in the following directory:
+
+.. code-block:: bash
+
+    raffle/example/python_pkg/graphite
+
+First, we need to establish an initial database of structures.
+This will be used to initialise the generalised distribution functions, which will inform the placement of atoms in the generated structures.
+Here, we detail two routes: 1) using an ASE object of bulk graphite and 2) using the Materials Project database.
+
+Here, we will simply use the ASE object of bulk graphite.
+If you wish to use the Materials Project database, please refer to the `Databases tutorial <databases>`_.
+
+First, we must import the required packages:
+
+.. code-block:: python
+
+    from ase import Atoms
+    from raffle.generator import raffle_generator
+    from mace.calculators import mace_mp
+    import numpy as np
+
+Next, we need to set up the RAFFLE generator and the calculator to calculate the energies of the structures.
+In this example, we use the MACE-MP-0 calculator:
+
+.. code-block:: python
+
+    generator = raffle_generator()
+
+    calc = mace_mp(model="medium", dispersion=False, default_dtype="float32", device='cpu')
+
+
+Then, we need to create the host structure.
+This is the base structure that will be added to in order to generate the structures.
+
+.. code-block:: python
+
+    host = Atoms(
+            "C6",
+            scaled_positions=[
+                [0.000000000, 0.000000000, 0.125000000],
+                [0.000000000, 0.000000000, 0.375000000],
+                [0.000000000, 0.000000000, 0.875000000],
+                [0.333333333, 0.666666667, 0.125000000],
+                [0.666666667, 0.333333333, 0.375000000],
+                [0.666666667, 0.333333333, 0.875000000],
+            ], cell=[
+                2.4672912616, 2.4672912616, 15.606146000
+            ], pbc=True
+        )
+
+    host.calc = calc
+    generator.set_host(host)
+
+Next, we need to set any variables we want to use for the distributions.
+None should be required other than the element energies, which will be used as reference energies to calculate the formation energies of the structures.
+However, for completeness, we will set other variables here.
+Because there is only one element in this structure, the reference energy is meaningless, so we set it to ``0.0`` here.
+This is due to all structures will have the same stoichiometry and thus use the exact same reference energy, so choice of this value is a global shift.
+
+.. code-block:: python
+
+    generator.distributions.set_element_energies( { 'C': 0.0 } )
+    generator.distributions.set_kBT(0.2)
+    generator.distributions.set_width([0.025, np.pi/200.0, np.pi/200.0])
+    generator.distributions.set_radius_distance_tol([1.5, 2.5, 3.0, 6.0])
+
+Now we need to generate the database of structures.
+We will provide bulk diamond as the only database entry here.
+
+.. code-block:: python
+
+    database = []
+    database.append(
+        Atoms("C4", 
+            positions=[
+                [0.0, 0.0, 1.95076825],
+                [0.0, 0.0, 5.85230475],
+                [1.2336456308015413, 0.7122456370278755, 1.95076825],
+                [1.2336456308015415, -0.7122456370278757, 5.85230475]
+            ], cell=[
+                [1.2336456308015413, -2.1367369110836267, 0.0], 
+                [1.2336456308015413,  2.1367369110836267, 0.0],
+                [0.0, 0.0, 7.803073]
+            ], pbc=True
+        )
+    )
+
+This database will now be used to initialise the generalised distribution functions in RAFFLE.
+
+.. code-block:: python
+
+    generator.distributions.create(database)
+
+Finally, we can set the grid on which atom searches are performed (this grid is applied to the host cell).
+By default, the grid is generated using a spacing of 0.1 Å.
+
+.. code-block:: python
+
+    generator.set_grid(grid_spacing=0.1, grid_offset=[0.0, 0.0, 0.0])
+
+We are now ready to generate structures using the database of structures.
+
+.. code-block:: python
+
+    num_structures_old = 0
+    structures, status = generator.generate(
+        num_structures = 1,
+        stoichiometry = { 'C': 2 },
+        method_ratio = {"void":0.0001, "min":1.0},
+        calc = calc
+    )
+
+We should now have a structure of layered graphite.
+This structure can be visualised using the ASE package.
+But this can also be verified energetically.
+The generated structure should have double the energy of bulk graphite, found in ```database[0]```.
\ No newline at end of file
diff --git a/docs/source/tutorials/index.rst b/docs/source/tutorials/index.rst
index e0bb1311..6f458851 100644
--- a/docs/source/tutorials/index.rst
+++ b/docs/source/tutorials/index.rst
@@ -18,6 +18,5 @@ Whilst RAFFLE is a random sturcture search package designed primarlily for inter
    databases_tutorial
    diamond_tutorial
    graphite_tutorial
-   BaTiO3_tutorial
-   perovskites_tutorial
-   C-MgO_tutorial
+   aluminium_tutorial
+   Si-Ge_tutorial
diff --git a/docs/source/tutorials/perovskites_tutorial.rst b/docs/source/tutorials/perovskites_tutorial.rst
deleted file mode 100644
index d378b21f..00000000
--- a/docs/source/tutorials/perovskites_tutorial.rst
+++ /dev/null
@@ -1,5 +0,0 @@
-.. perovskites:
-
-====================================================
-BaTiO\ :sub:`3`\|SrTiO\ :sub:`3` interface tutorial
-====================================================
\ No newline at end of file
diff --git a/example/python_pkg/Al_learn/learn.py b/example/python_pkg/Al_learn/learn.py
index e2010ff3..cc793de5 100644
--- a/example/python_pkg/Al_learn/learn.py
+++ b/example/python_pkg/Al_learn/learn.py
@@ -63,9 +63,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct
 
     # set up the hosts
     crystal_structures = [
-        'orthorhombic', 'diamond',
-        'bct', 'sc',
-        'fcc', 'bcc', 'hcp',
+        'orthorhombic', 'hcp',
     ]
     hosts = []
     for crystal_structure in crystal_structures:
diff --git a/example/python_pkg/Si-Ge_learn/learn.py b/example/python_pkg/Si-Ge_learn/learn.py
index b7cbb658..4697fdd8 100644
--- a/example/python_pkg/Si-Ge_learn/learn.py
+++ b/example/python_pkg/Si-Ge_learn/learn.py
@@ -1,4 +1,4 @@
-from chgnet.model import CHGNetCalculator
+from mace.calculators import mace_mp
 from raffle.generator import raffle_generator
 from ase import build
 from ase.optimize import FIRE
@@ -57,8 +57,8 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct
 if __name__ == "__main__":
 
     # set up the calculator
-    calc_params = {}
-    calc = CHGNetCalculator()
+    calc_params = { 'model': "../mace-mpa-0-medium.model" }
+    calc = mace_mp(**calc_params)
 
     # set up the hosts
     hosts = []
@@ -136,7 +136,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct
     num_structures_old = 0
     optimise_structure = True
     # start the iterations, loop over the hosts, then repeat the process X times
-    for ival in range(20):
+    for ival in range(40):
         for host in hosts:
             print("setting host")
             generator.set_host(host)
@@ -150,7 +150,7 @@ def process_structure(i, atoms, num_structures_old, calc_params, optimise_struct
                 num_structures = 5,
                 stoichiometry = { 'Si': 16, 'Ge': 16 },
                 seed = seed*1000+iter,
-                method_ratio = {"void": 0.3, "rand": 0.01, "walk": 0.3, "grow": 0.0, "min": 1.0},
+                method_ratio = {"void": 0.1, "rand": 0.01, "walk": 0.25, "grow": 0.25, "min": 1.0},
                 verbose = 0,
             )
 
diff --git a/example/python_pkg/diamond/POSCAR_host b/example/python_pkg/diamond/POSCAR_host
new file mode 100644
index 00000000..11242071
--- /dev/null
+++ b/example/python_pkg/diamond/POSCAR_host
@@ -0,0 +1,16 @@
+C8
+    1.000000000
+    3.560745109    0.000000000    0.000000000
+    0.000000000    3.560745109    0.000000000
+    0.000000000    0.000000000    7.121490218
+C
+8
+Direct
+    0.000000000    0.000000000    0.000000000
+    0.500000000    0.500000000    0.000000000
+    0.500000000    0.000000000    0.250000000
+    0.000000000    0.500000000    0.250000000
+    0.250000000    0.250000000    0.125000000
+    0.750000000    0.750000000    0.125000000
+    0.750000000    0.250000000    0.375000000
+    0.250000000    0.750000000    0.375000000
diff --git a/example/python_pkg/diamond/POSCAR_host_3x3 b/example/python_pkg/diamond/POSCAR_host_3x3
new file mode 100644
index 00000000..72b5deab
--- /dev/null
+++ b/example/python_pkg/diamond/POSCAR_host_3x3
@@ -0,0 +1,80 @@
+C8
+    1.000000000
+   14.242980436    0.000000000    0.000000000
+    0.000000000   14.242980436    0.000000000
+    0.000000000    0.000000000    7.121490218
+C    
+72
+Direct
+    0.000000000    0.000000000    0.000000000
+    0.333333333    0.000000000    0.000000000
+    0.666666667    0.000000000    0.000000000
+    0.000000000    0.333333333    0.000000000
+    0.333333333    0.333333333    0.000000000
+    0.666666667    0.333333333    0.000000000
+    0.000000000    0.666666667    0.000000000
+    0.333333333    0.666666667    0.000000000
+    0.666666667    0.666666667    0.000000000
+    0.166666667    0.166666667    0.000000000
+    0.500000000    0.166666667    0.000000000
+    0.833333333    0.166666667    0.000000000
+    0.166666667    0.500000000    0.000000000
+    0.500000000    0.500000000    0.000000000
+    0.833333333    0.500000000    0.000000000
+    0.166666667    0.833333333    0.000000000
+    0.500000000    0.833333333    0.000000000
+    0.833333333    0.833333333    0.000000000
+    0.166666667    0.000000000    0.250000000
+    0.500000000    0.000000000    0.250000000
+    0.833333333    0.000000000    0.250000000
+    0.166666667    0.333333333    0.250000000
+    0.500000000    0.333333333    0.250000000
+    0.833333333    0.333333333    0.250000000
+    0.166666667    0.666666667    0.250000000
+    0.500000000    0.666666667    0.250000000
+    0.833333333    0.666666667    0.250000000
+    0.000000000    0.166666667    0.250000000
+    0.333333333    0.166666667    0.250000000
+    0.666666667    0.166666667    0.250000000
+    0.000000000    0.500000000    0.250000000
+    0.333333333    0.500000000    0.250000000
+    0.666666667    0.500000000    0.250000000
+    0.000000000    0.833333333    0.250000000
+    0.333333333    0.833333333    0.250000000
+    0.666666667    0.833333333    0.250000000
+    0.083333333    0.083333333    0.125000000
+    0.416666667    0.083333333    0.125000000
+    0.750000000    0.083333333    0.125000000
+    0.083333333    0.416666667    0.125000000
+    0.416666667    0.416666667    0.125000000
+    0.750000000    0.416666667    0.125000000
+    0.083333333    0.750000000    0.125000000
+    0.416666667    0.750000000    0.125000000
+    0.750000000    0.750000000    0.125000000
+    0.250000000    0.250000000    0.125000000
+    0.583333333    0.250000000    0.125000000
+    0.916666667    0.250000000    0.125000000
+    0.250000000    0.583333333    0.125000000
+    0.583333333    0.583333333    0.125000000
+    0.916666667    0.583333333    0.125000000
+    0.250000000    0.916666667    0.125000000
+    0.583333333    0.916666667    0.125000000
+    0.916666667    0.916666667    0.125000000
+    0.250000000    0.083333333    0.375000000
+    0.583333333    0.083333333    0.375000000
+    0.916666667    0.083333333    0.375000000
+    0.250000000    0.416666667    0.375000000
+    0.583333333    0.416666667    0.375000000
+    0.916666667    0.416666667    0.375000000
+    0.250000000    0.750000000    0.375000000
+    0.583333333    0.750000000    0.375000000
+    0.916666667    0.750000000    0.375000000
+    0.083333333    0.250000000    0.375000000
+    0.416666667    0.250000000    0.375000000
+    0.750000000    0.250000000    0.375000000
+    0.083333333    0.583333333    0.375000000
+    0.416666667    0.583333333    0.375000000
+    0.750000000    0.583333333    0.375000000
+    0.083333333    0.916666667    0.375000000
+    0.416666667    0.916666667    0.375000000
+    0.750000000    0.916666667    0.375000000
diff --git a/example/python_pkg/graphite/POSCAR_host_graphene b/example/python_pkg/graphite/POSCAR_host_graphene
new file mode 100644
index 00000000..f5d05414
--- /dev/null
+++ b/example/python_pkg/graphite/POSCAR_host_graphene
@@ -0,0 +1,26 @@
+C4
+    1.000000000
+    3.700936892   -6.410210733    0.000000000
+    3.700936892    6.410210733    0.000000000
+    0.000000000    0.000000000    7.803073000
+C    
+18
+Direct
+    0.000000000    0.000000000    0.250000000
+    0.333333333    0.000000000    0.250000000
+    0.666666667    0.000000000    0.250000000
+    0.000000000    0.333333333    0.250000000
+    0.333333333    0.333333333    0.250000000
+    0.666666667    0.333333333    0.250000000
+    0.000000000    0.666666667    0.250000000
+    0.333333333    0.666666667    0.250000000
+    0.666666667    0.666666667    0.250000000
+    0.111111111    0.222222222    0.250000000
+    0.444444444    0.222222222    0.250000000
+    0.777777778    0.222222222    0.250000000
+    0.111111111    0.555555556    0.250000000
+    0.444444444    0.555555556    0.250000000
+    0.777777778    0.555555556    0.250000000
+    0.111111111    0.888888889    0.250000000
+    0.444444444    0.888888889    0.250000000
+    0.777777778    0.888888889    0.250000000
\ No newline at end of file
diff --git a/example/python_pkg/graphite/POSCAR_host_graphite_vacancy b/example/python_pkg/graphite/POSCAR_host_graphite_vacancy
new file mode 100644
index 00000000..4fb28784
--- /dev/null
+++ b/example/python_pkg/graphite/POSCAR_host_graphite_vacancy
@@ -0,0 +1,44 @@
+C4
+    1.000000000
+    3.700936892   -6.410210733    0.000000000
+    3.700936892    6.410210733    0.000000000
+    0.000000000    0.000000000    7.803073000
+C    
+35
+Direct
+    0.000000000    0.000000000    0.250000000
+    0.333333333    0.000000000    0.250000000
+    0.666666667    0.000000000    0.250000000
+    0.000000000    0.333333333    0.250000000
+    0.333333333    0.333333333    0.250000000
+    0.666666667    0.333333333    0.250000000
+    0.000000000    0.666666667    0.250000000
+    0.333333333    0.666666667    0.250000000
+    0.666666667    0.666666667    0.250000000
+    0.000000000    0.000000000    0.750000000
+    0.333333333    0.000000000    0.750000000
+    0.666666667    0.000000000    0.750000000
+    0.000000000    0.333333333    0.750000000
+    0.333333333    0.333333333    0.750000000
+    0.666666667    0.333333333    0.750000000
+    0.000000000    0.666666667    0.750000000
+    0.333333333    0.666666667    0.750000000
+    0.666666667    0.666666667    0.750000000
+    0.111111111    0.222222222    0.250000000
+    0.444444444    0.222222222    0.250000000
+    0.777777778    0.222222222    0.250000000
+    0.111111111    0.555555556    0.250000000
+    0.444444444    0.555555556    0.250000000
+    0.777777778    0.555555556    0.250000000
+    0.111111111    0.888888889    0.250000000
+    0.444444444    0.888888889    0.250000000
+    0.777777778    0.888888889    0.250000000
+    0.222222222    0.111111111    0.750000000
+    0.555555556    0.111111111    0.750000000
+    0.888888889    0.111111111    0.750000000
+    0.222222222    0.444444444    0.750000000
+    0.555555556    0.444444444    0.750000000
+    0.888888889    0.444444444    0.750000000
+    0.222222222    0.777777778    0.750000000
+    0.555555556    0.777777778    0.750000000
+
diff --git a/example/python_pkg/graphite/POSCAR_host_missing_layer b/example/python_pkg/graphite/POSCAR_host_missing_layer
new file mode 100644
index 00000000..024e65f9
--- /dev/null
+++ b/example/python_pkg/graphite/POSCAR_host_missing_layer
@@ -0,0 +1,15 @@
+C
+    1.000000000
+    1.233645631   -2.136736911    0.000000000
+    1.233645631    2.136736911    0.000000000
+    0.000000000    0.000000000   15.606146000
+C    
+6
+Direct
+    0.000000000    0.000000000    0.125000000
+    0.000000000    0.000000000    0.375000000
+    0.000000000    0.000000000    0.875000000
+    0.333333333    0.666666667    0.125000000
+    0.666666667    0.333333333    0.375000000
+    0.666666667    0.333333333    0.875000000
+
diff --git a/src/fortran/lib/mod_dist_calcs.f90 b/src/fortran/lib/mod_dist_calcs.f90
index c92621f9..9bade3df 100644
--- a/src/fortran/lib/mod_dist_calcs.f90
+++ b/src/fortran/lib/mod_dist_calcs.f90
@@ -84,8 +84,8 @@ function get_min_dist(basis,loc,lignore_close,axis,labove,lreal,tol, &
     do js = 1, basis%nspec
        atmloop: do ja = 1, basis%spec(js)%num
           if(present(ignore_list))then
-             do i = 1, size(ignore_list,1), 1
-                if(all(ignore_list(i,:).eq.[js,ja])) cycle atmloop
+             do i = 1, size(ignore_list,2), 1
+                if(all(ignore_list(:,i).eq.[js,ja])) cycle atmloop
              end do
           end if
           vdtmp1 = basis%spec(js)%atom(ja,:3) - loc
@@ -181,8 +181,8 @@ pure function get_min_dist_between_point_and_species( &
     dist = huge(0._real32)
     atom_loop: do ia = 1,basis%spec(species)%num
        if(present(ignore_list))then
-          do i = 1, size(ignore_list,1), 1
-             if(all(ignore_list(i,:).eq.[species,ia])) cycle atom_loop
+          do i = 1, size(ignore_list,2), 1
+             if(all(ignore_list(:,i).eq.[species,ia])) cycle atom_loop
           end do
        end if
        vec = loc - basis%spec(species)%atom(ia,:3)
diff --git a/src/fortran/lib/mod_evaluator.f90 b/src/fortran/lib/mod_evaluator.f90
index 6956b2e0..cb786b5c 100644
--- a/src/fortran/lib/mod_evaluator.f90
+++ b/src/fortran/lib/mod_evaluator.f90
@@ -107,8 +107,8 @@ function evaluate_point( distribs_container, &
        atom_loop: do ia = 1, basis%spec(is)%num
           ! Check if the atom is in the ignore list
           ! If it is, skip the atom.
-          do i = 1, size(atom_ignore_list,dim=1), 1
-             if(all(atom_ignore_list(i,:).eq.[is,ia])) cycle atom_loop
+          do i = 1, size(atom_ignore_list,dim=2), 1
+             if(all(atom_ignore_list(:,i).eq.[is,ia])) cycle atom_loop
           end do
           associate( position_store => [ basis%spec(is)%atom(ia,1:3) ] )
              bondlength = modu( matmul(position - position_store, basis%lat) )
diff --git a/src/fortran/lib/mod_generator.f90 b/src/fortran/lib/mod_generator.f90
index ea862d0d..35f59968 100644
--- a/src/fortran/lib/mod_generator.f90
+++ b/src/fortran/lib/mod_generator.f90
@@ -298,7 +298,7 @@ end subroutine reset_bounds
 !###############################################################################
   subroutine generate(this, num_structures, &
        stoichiometry, method_ratio, seed, settings_out_file, &
-       verbose &
+       verbose, exit_code &
   )
     !! Generate random structures.
     !!
@@ -323,10 +323,14 @@ subroutine generate(this, num_structures, &
     !! File to print the settings to.
     integer, intent(in), optional :: verbose
     !! Verbosity level.
+    integer, intent(out), optional :: exit_code
+    !! Exit code.
 
     ! Local variables
     integer :: i, j, k, istructure, num_structures_old, num_structures_new
     !! Loop counters.
+    integer :: exit_code_
+    !! Exit code.
     integer :: num_seed
     !! Number of seeds for the random number generator.
     integer :: num_insert_atoms, num_insert_species
@@ -356,6 +360,7 @@ subroutine generate(this, num_structures, &
     !---------------------------------------------------------------------------
     ! Set the verbosity level
     !---------------------------------------------------------------------------
+    exit_code_ = 0
     verbose_ = 0
     if(present(verbose)) verbose_ = verbose
     if(verbose_ .eq. 0) suppress_warnings = .true.
@@ -467,7 +472,7 @@ subroutine generate(this, num_structures, &
     ! ... basis_template
     !---------------------------------------------------------------------------
     if(verbose_.gt.0) write(*,*) "Generating placement list"
-    allocate(placement_list(num_insert_atoms,2))
+    allocate(placement_list(2, num_insert_atoms))
     k = 0
     spec_loop1: do i = 1, basis_template%nspec
        success = .false.
@@ -481,15 +486,15 @@ subroutine generate(this, num_structures, &
        if(i.gt.this%host%nspec)then
           do j = 1, basis_template%spec(i)%num
              k = k + 1
-             placement_list(k,1) = i
-             placement_list(k,2) = j
+             placement_list(1,k) = i
+             placement_list(2,k) = j
           end do
        else
           do j = 1, basis_template%spec(i)%num
              if(j.le.this%host%spec(i)%num) cycle
              k = k + 1
-             placement_list(k,1) = i
-             placement_list(k,2) = j
+             placement_list(1,k) = i
+             placement_list(2,k) = j
           end do
        end if
     end do spec_loop1
@@ -509,16 +514,24 @@ subroutine generate(this, num_structures, &
                  basis_template, &
                  placement_list, &
                  method_rand_limit, &
-                 verbose_ &
+                 verbose_, &
+                 exit_code_ &
             ) &
        )
        this%num_structures = istructure
 
     end do structure_loop
-    if(verbose_ .gt. 0) write(*,*) "Finished generating structures"
+    if(verbose_ .gt. 0 .and. exit_code_ .eq. 0) &
+         write(*,*) "Finished generating structures"
 
     if(verbose_ .eq. 0) suppress_warnings = .true.
 
+    if(present(exit_code))then
+       exit_code = exit_code_
+    elseif(exit_code_ .ne. 0)then
+       call stop_program("Error generating structures", exit_code_)
+    end if
+
   end subroutine generate
 !###############################################################################
 
@@ -527,7 +540,8 @@ end subroutine generate
   function generate_structure( &
        this, &
        basis_initial, &
-       placement_list, method_rand_limit, verbose &
+       placement_list, method_rand_limit, verbose, &
+       exit_code &
   ) result(basis)
     !! Generate a single random structure.
     !!
@@ -551,6 +565,8 @@ function generate_structure( &
     !! Generated basis.
     integer, intent(in) :: verbose
     !! Verbosity level.
+    integer, intent(inout) :: exit_code
+    !! Exit code.
 
     ! Local variables
     integer :: iplaced, void_ticker
@@ -561,6 +577,8 @@ function generate_structure( &
     !! Random number.
     logical :: viable
     !! Boolean for viable placement.
+    logical :: placement_aborted
+    !! Boolean for aborted placement.
     integer, dimension(size(placement_list,1),size(placement_list,2)) :: &
          placement_list_shuffled
     !! Shuffled placement list.
@@ -593,13 +611,13 @@ function generate_structure( &
     ! shuffle the placement list
     !---------------------------------------------------------------------------
     placement_list_shuffled = placement_list
-    call shuffle(placement_list_shuffled,1)
+    call shuffle(placement_list_shuffled,2)
 
 
     !---------------------------------------------------------------------------
     ! generate species index list to add
     !---------------------------------------------------------------------------
-    species_index_list = placement_list_shuffled(:,1)
+    species_index_list = placement_list_shuffled(1,:)
     call set(species_index_list)
 
 
@@ -627,6 +645,7 @@ function generate_structure( &
     iplaced = 0
     void_ticker = 0
     viable = .false.
+    placement_aborted = .false.
     placement_loop: do while (iplaced.lt.num_insert_atoms)
        !------------------------------------------------------------------------
        ! check if there are any viable gridpoints remaining
@@ -638,19 +657,15 @@ function generate_structure( &
                     this%distributions, &
                     basis, &
                     species_index_list, &
-                    [ placement_list_shuffled(iplaced,:) ], &
+                    [ placement_list_shuffled(:,iplaced) ], &
                     [ this%distributions%bond_info(:)%radius_covalent ], &
-                    placement_list_shuffled(iplaced+1:,:) &
+                    placement_list_shuffled(:,iplaced+1:) &
                )
           if(.not.allocated(gridpoint_viability))then
              if(abs(method_rand_limit(4)).lt.1.E-6)then
-                call stop_program( &
-                     "No viable gridpoints" // &
-                     achar(13) // achar(10) // &
-                     "No placement methods available" &
-                )
-                return
-             else if( &
+                placement_aborted = .true.
+                exit placement_loop
+             elseif( &
                   abs( &
                        method_rand_limit(5) - method_rand_limit(4) &
                   ) .gt. 1.E-6 &
@@ -673,45 +688,42 @@ function generate_structure( &
        call random_number(rtmp1)
        if(rtmp1.le.method_rand_limit(1)) then
           if(verbose.gt.0) write(*,*) "Add Atom Void"
-          point = place_method_void( this%grid, &
-               this%grid_offset, &
-               this%bounds, &
+          point = place_method_void( &
+               gridpoint_viability, &
                basis, &
-               placement_list_shuffled(iplaced+1:,:), viable &
+               placement_list_shuffled(:,iplaced+1:), viable &
           )
-       else if(rtmp1.le.method_rand_limit(2)) then
+       elseif(rtmp1.le.method_rand_limit(2)) then
           if(verbose.gt.0) write(*,*) "Add Atom Random"
           point = place_method_rand( &
                this%distributions, &
                this%bounds, &
                basis, &
-               placement_list_shuffled(iplaced+1:,:), &
+               placement_list_shuffled(:,iplaced+1:), &
                [ this%distributions%bond_info(:)%radius_covalent ], &
                this%max_attempts, &
                viable &
           )
-          if(.not. viable) cycle placement_loop
-       else if(rtmp1.le.method_rand_limit(3)) then
+       elseif(rtmp1.le.method_rand_limit(3)) then
           if(verbose.gt.0) write(*,*) "Add Atom Walk"
           point = place_method_walk( &
                this%distributions, &
                this%bounds, &
                basis, &
-               placement_list_shuffled(iplaced+1:,:), &
+               placement_list_shuffled(:,iplaced+1:), &
                [ this%distributions%bond_info(:)%radius_covalent ], &
                this%max_attempts, &
                this%walk_step_size_coarse, this%walk_step_size_fine, &
                viable &
           )
-          if(.not. viable) void_ticker = void_ticker + 1
-       else if(rtmp1.le.method_rand_limit(4)) then
+       elseif(rtmp1.le.method_rand_limit(4)) then
           if(iplaced.eq.0)then
              if(verbose.gt.0) write(*,*) "Add Atom Random (growth seed)"
              point = place_method_rand( &
                   this%distributions, &
                   this%bounds, &
                   basis, &
-                  placement_list_shuffled(iplaced+1:,:), &
+                  placement_list_shuffled(:,iplaced+1:), &
                   [ this%distributions%bond_info(:)%radius_covalent ], &
                   this%max_attempts, &
                   viable &
@@ -720,38 +732,29 @@ function generate_structure( &
              if(verbose.gt.0) write(*,*) "Add Atom Growth"
              point = place_method_growth( &
                   this%distributions, &
-                  basis%spec(placement_list_shuffled(iplaced,1))%atom( &
-                       placement_list_shuffled(iplaced,2),:3 &
+                  basis%spec(placement_list_shuffled(1,iplaced))%atom( &
+                       placement_list_shuffled(2,iplaced),:3 &
                   ), &
-                  placement_list_shuffled(iplaced,1), &
+                  placement_list_shuffled(1,iplaced), &
                   this%bounds, &
                   basis, &
-                  placement_list_shuffled(iplaced+1:,:), &
+                  placement_list_shuffled(:,iplaced+1:), &
                   [ this%distributions%bond_info(:)%radius_covalent ], &
                   this%max_attempts, &
                   this%walk_step_size_coarse, this%walk_step_size_fine, &
                   viable &
              )
           end if
-          if(.not. viable) void_ticker = void_ticker + 1
-       else if(rtmp1.le.method_rand_limit(5)) then
+       elseif(rtmp1.le.method_rand_limit(5)) then
           if(verbose.gt.0) write(*,*) "Add Atom Minimum"
           point = place_method_min( gridpoint_viability, &
-               placement_list_shuffled(iplaced+1,1), &
+               placement_list_shuffled(1,iplaced+1), &
                species_index_list, &
                viable &
           )
           if(.not. viable .and. abs(method_rand_limit(4)).lt.1.E-6)then
-             write(stop_msg,*) &
-                  "No viable gridpoints" // &
-                  achar(13) // achar(10) // &
-                  "Min method is the only method, but cannot place another &
-                  &atom" // &
-                  achar(13) // achar(10) // &
-                  "Species to place now: ", &
-                  basis%spec(placement_list_shuffled(iplaced+1,1))%name
-             call stop_program(stop_msg)
-             return
+             placement_aborted = .true.
+             exit placement_loop
           elseif(.not. viable)then
              deallocate(gridpoint_viability)
              write(warn_msg, '("No more viable gridpoints")')
@@ -767,25 +770,30 @@ function generate_structure( &
        ! check if the placement method returned a viable point
        ! if not, cycle the loop
        !------------------------------------------------------------------------
-       if(.not. viable) then
-          if(void_ticker.gt.10) &
-               point = place_method_void( &
-                    this%grid, this%grid_offset, this%bounds, basis, &
-                    placement_list_shuffled(iplaced+1:,:), viable &
-               )
-          void_ticker = 0
+       if(.not. viable)then
+          void_ticker = void_ticker + 1
+          if(void_ticker.gt.10.and..not.allocated(gridpoint_viability))then
+             placement_aborted = .true.
+             exit placement_loop
+          elseif(void_ticker.gt.10)then
+             point = place_method_void( &
+                  gridpoint_viability, basis, &
+                  placement_list_shuffled(:,iplaced+1:), viable &
+             )
+             void_ticker = 0
+          end if
           if(.not.viable) cycle placement_loop
        end if
        !------------------------------------------------------------------------
        ! place the atom and update the image atoms in the basis
        !------------------------------------------------------------------------
        iplaced = iplaced + 1
-       basis%spec(placement_list_shuffled(iplaced,1))%atom( &
-            placement_list_shuffled(iplaced,2),:3) = point(:3)
+       basis%spec(placement_list_shuffled(1,iplaced))%atom( &
+            placement_list_shuffled(2,iplaced),:3) = point(:3)
        call basis%update_images( &
             max_bondlength = this%distributions%cutoff_max(1), &
-            is = placement_list_shuffled(iplaced,1), &
-            ia = placement_list_shuffled(iplaced,2) &
+            is = placement_list_shuffled(1,iplaced), &
+            ia = placement_list_shuffled(2,iplaced) &
        )
        if(verbose.gt.0)then
           write(*,'(A)',ADVANCE='NO') achar(13)
@@ -793,6 +801,16 @@ function generate_structure( &
        end if
 
     end do placement_loop
+    
+    if(placement_aborted)then
+       call stop_program( &
+            "Placement routine aborted, not all atoms placed", &
+            block_stop = .true. &
+       )
+       exit_code = 1
+       call basis%remove_atoms(placement_list_shuffled(:,iplaced+1:))
+    end if
+
     if(allocated(gridpoint_viability)) deallocate(gridpoint_viability)
 
   end function generate_structure
@@ -922,7 +940,7 @@ function evaluate(this, basis) result(viability)
     ! Local variables
     integer :: is, ia, species
     !! Loop indices.
-    integer, dimension(1,2) :: atom_ignore
+    integer, dimension(2,1) :: atom_ignore
     !! Atom to ignore.
     type(extended_basis_type) :: basis_extd
     !! Extended basis for the structure to evaluate.
@@ -949,7 +967,7 @@ function evaluate(this, basis) result(viability)
           return
        end if
        do ia = 1, basis%spec(is)%num
-          atom_ignore(1,:) = [is,ia]
+          atom_ignore(:,1) = [is,ia]
           viability = viability + &
                evaluate_point( this%distributions, &
                     [ basis%spec(is)%atom(ia,1:3) ], &
diff --git a/src/fortran/lib/mod_geom_extd.f90 b/src/fortran/lib/mod_geom_extd.f90
index 3c27f13d..3b1d4dd0 100644
--- a/src/fortran/lib/mod_geom_extd.f90
+++ b/src/fortran/lib/mod_geom_extd.f90
@@ -92,8 +92,8 @@ subroutine create_images(this, max_bondlength, atom_ignore_list)
        image_species(is)%radius = this%spec(is)%radius
        image_species(is)%name = this%spec(is)%name
        atom_loop: do ia = 1, this%spec(is)%num
-          do i = 1, size(atom_ignore_list_,1), 1
-             if(all(atom_ignore_list_(i,:).eq.[is,ia])) cycle atom_loop
+          do i = 1, size(atom_ignore_list_,2), 1
+             if(all(atom_ignore_list_(:,i).eq.[is,ia])) cycle atom_loop
           end do
           do i=-amax,amax+1,1
              vtmp1(1) = this%spec(is)%atom(ia,1) + real(i, real32)
diff --git a/src/fortran/lib/mod_geom_rw.f90 b/src/fortran/lib/mod_geom_rw.f90
index 720630d9..f4eecf3f 100644
--- a/src/fortran/lib/mod_geom_rw.f90
+++ b/src/fortran/lib/mod_geom_rw.f90
@@ -71,6 +71,10 @@ module raffle__geom_rw
      !! Procedure to copy the basis.
      procedure, pass(this) :: get_lattice_constants
      !! Procedure to get the lattice constants of the basis.
+     procedure, pass(this) :: remove_atom
+     !! Procedure to remove an atom from the basis.
+     procedure, pass(this) :: remove_atoms
+     !! Procedure to remove atoms from the basis.
   end type basis_type
 
 
@@ -1406,6 +1410,137 @@ end subroutine copy
 !###############################################################################
 
 
+!###############################################################################
+  subroutine remove_atom(this, ispec, iatom)
+    !! Remove an atom from the basis.
+    implicit none
+
+    ! Arguments
+    class(basis_type), intent(inout) :: this
+    !! Parent. The basis.
+    integer, intent(in) :: ispec, iatom
+    !! The species and atom to remove.
+
+    ! Local variables
+    integer :: i
+    !! Loop index.
+    real(real32), dimension(:,:), allocatable :: atom
+    !! Temporary array to store the atomic positions.
+
+
+    !---------------------------------------------------------------------------
+    ! remove atom from basis
+    !---------------------------------------------------------------------------
+    do i = 1, this%nspec
+       if(i.eq.ispec)then
+          if(iatom.gt.this%spec(i)%num)then
+             call stop_program("Atom to remove does not exist")
+             return
+          end if
+          allocate(atom(this%spec(i)%num-1,size(this%spec(i)%atom,2)))
+          atom(1:iatom-1:1,:) = this%spec(i)%atom(1:iatom-1:1,:)
+          atom(iatom:this%spec(i)%num-1:1,:) = &
+               this%spec(i)%atom(iatom+1:this%spec(i)%num:1,:)
+          this%spec(i)%atom = atom
+          deallocate(atom)
+          this%spec(i)%num = this%spec(i)%num - 1
+          this%natom = this%natom - 1
+          if(this%spec(i)%num.eq.0)then
+             deallocate(this%spec(i)%atom)
+             if(this%nspec.eq.0)then
+                deallocate(this%spec)
+                this%lcart = .true.
+                this%sysname = ""
+                this%energy = 0._real32
+                this%lat = 0._real32
+                this%pbc = .true.
+             end if
+          end if
+       end if
+    end do
+
+  end subroutine remove_atom
+!###############################################################################
+
+
+!###############################################################################
+  subroutine remove_atoms(this, atoms)
+    !! Remove atoms from the basis.
+    use raffle__misc, only: swap
+    implicit none
+
+    ! Arguments
+    class(basis_type), intent(inout) :: this
+    !! Parent. The basis.
+    integer, dimension(:,:), intent(in) :: atoms
+    !! The atoms to remove (2, number of atoms to remove)
+    !! 1st value of 1st dimension is the species number
+    !! 2nd value of 1st dimension is the atom number
+    !! 2nd dimension is the number of atoms to remove
+
+    ! Local variables
+    integer :: is, ia, i
+    !! Loop index.
+    integer :: n, m, start_idx, end_idx, loc
+    !! Index variables.
+    integer :: num_species
+    !! The number of species.
+    integer, dimension(:,:), allocatable :: atoms_ordered
+    !! The atoms to remove ordered by species and atom
+    real(real32), dimension(:,:), allocatable :: atom
+    !! Temporary array to store the atomic positions.
+
+
+    !---------------------------------------------------------------------------
+    ! reorder atoms to remove
+    !---------------------------------------------------------------------------
+    allocate(atoms_ordered, source=atoms)
+    n = size(atoms_ordered, 1)
+    m = size(atoms_ordered, 2)
+
+    do i = 1, m
+       loc = maxloc(atoms_ordered(1, i:n), dim=1) + i - 1
+       if (loc .ne. i) then
+          call swap(atoms_ordered(1, i), atoms_ordered(1, loc))
+          call swap(atoms_ordered(2, i), atoms_ordered(2, loc))
+       end if
+    end do
+    num_species = this%nspec
+    do is = 1, num_species
+       start_idx = findloc(atoms_ordered(1, :), is, dim=1)
+       end_idx   = findloc(atoms_ordered(1, :), is, dim=1, back=.true.)
+       if (start_idx .eq. 0) cycle
+       do ia = start_idx, end_idx, 1
+          loc = maxloc( &
+               atoms_ordered(2, ia:end_idx), &
+               dim=1 &
+          ) + ia - 1
+          if (loc .ne. ia) then
+             call swap(atoms_ordered(1, ia), atoms_ordered(1, loc))
+             call swap(atoms_ordered(2, ia), atoms_ordered(2, loc))
+          end if
+       end do
+    end do
+
+
+    !---------------------------------------------------------------------------
+    ! remove atoms from basis
+    !---------------------------------------------------------------------------
+    do i = 1, size(atoms_ordered, 2)
+       call this%remove_atom(atoms_ordered(1, i), atoms_ordered(2, i))
+    end do
+
+    do is = 1, this%nspec
+       if (this%spec(is)%num .eq. 0) then
+          this%spec = [ this%spec(1:is-1), this%spec(is+1:) ]
+          this%nspec = this%nspec - 1
+       end if
+    end do
+
+  end subroutine remove_atoms
+!###############################################################################
+
+
 !###############################################################################
   subroutine get_element_properties(element, charge, mass, radius)
     !! Set the mass and charge of the element
diff --git a/src/fortran/lib/mod_io_utils.F90 b/src/fortran/lib/mod_io_utils.F90
index 9a7ee34a..d36e6d81 100644
--- a/src/fortran/lib/mod_io_utils.F90
+++ b/src/fortran/lib/mod_io_utils.F90
@@ -20,23 +20,32 @@ module raffle__io_utils
 contains
 
 !###############################################################################
-  subroutine stop_program(message, exit_code)
+  subroutine stop_program(message, exit_code, block_stop)
     !! Stop the program and print an error message.
     implicit none
     character(len=*), intent(in) :: message
     integer, intent(in), optional :: exit_code
+    logical, intent(in), optional :: block_stop
 
     integer :: exit_code_
+    logical :: block_stop_
 
     if(present(exit_code)) then
        exit_code_ = exit_code
     else
        exit_code_ = 1
     end if
+    if(present(block_stop)) then
+       block_stop_ = block_stop
+    else
+       block_stop_ = .false.
+    end if
 
     write(0,*) 'ERROR: ', trim(message)
-    if(.not.test_error_handling)then
-       stop exit_code_
+    if(.not.block_stop_)then
+       if(.not.test_error_handling) then
+          stop exit_code_
+       end if
     end if
   end subroutine stop_program
 !###############################################################################
diff --git a/src/fortran/lib/mod_misc.f90 b/src/fortran/lib/mod_misc.f90
index e0d6405b..83a3ec69 100644
--- a/src/fortran/lib/mod_misc.f90
+++ b/src/fortran/lib/mod_misc.f90
@@ -620,13 +620,14 @@ subroutine ishuffle(arr,dim,seed)
        iother = 2
        i2s=1;i2e=size(arr,dim=iother)
        j2s=1;j2e=size(arr,dim=iother)
+       allocate(tlist(1,size(arr,dim=iother)))
     else
        iother = 1
        i1s=1;i1e=size(arr,dim=iother)
        j1s=1;j1e=size(arr,dim=iother)
+       allocate(tlist(size(arr,dim=iother),1))
     end if
     istart=1
-    allocate(tlist(1,size(arr,dim=iother)))
     do k = 1, 2
        do i = 1, n_data
           call random_number(r)
@@ -638,9 +639,9 @@ subroutine ishuffle(arr,dim,seed)
              i2s=i;i2e=i
              j2s=j;j2e=j
           end if
-          tlist(1:1,:) = arr(i1s:i1e,i2s:i2e)
+          tlist(:,:) = arr(i1s:i1e,i2s:i2e)
           arr(i1s:i1e,i2s:i2e) = arr(j1s:j1e,j2s:j2e)
-          arr(j1s:j1e,j2s:j2e) = tlist(1:1,:)
+          arr(j1s:j1e,j2s:j2e) = tlist(:,:)
        end do
     end do
 
diff --git a/src/fortran/lib/mod_place_methods.f90 b/src/fortran/lib/mod_place_methods.f90
index cb704429..bb3e17ab 100644
--- a/src/fortran/lib/mod_place_methods.f90
+++ b/src/fortran/lib/mod_place_methods.f90
@@ -32,7 +32,7 @@ module raffle__place_methods
 
 !###############################################################################
   function place_method_void( &
-       grid, grid_offset, bounds, basis, atom_ignore_list, viable &
+       points, basis, atom_ignore_list, viable &
   ) result(point)
     !! VOID placement method.
     !!
@@ -41,14 +41,10 @@ function place_method_void( &
     implicit none
 
     ! Arguments
+    real(real32), dimension(:,:), intent(in) :: points
+    !! List of gridpoints to consider.
     type(extended_basis_type), intent(inout) :: basis
     !! Structure to add atom to.
-    integer, dimension(3), intent(in) :: grid
-    !! Number of gridpoints in each direction.
-    real(real32), dimension(3), intent(in) :: grid_offset
-    !! Offset for gridpoints.
-    real(real32), dimension(2,3), intent(in) :: bounds
-    !! Bounds of the unit cell.
     integer, dimension(:,:), intent(in) :: atom_ignore_list
     !! List of atoms to ignore (i.e. indices of atoms not yet placed).
     logical, intent(out) :: viable
@@ -75,21 +71,14 @@ function place_method_void( &
     !---------------------------------------------------------------------------
     best_location = 0._real32
     best_location_bond = -huge(1._real32)
-    do i = 0, grid(1) - 1, 1
-       do j = 0, grid(2) - 1, 1
-          do k = 0, grid(3) - 1, 1
-             tmpvector = bounds(1,:) + &
-                  ( bounds(2,:) - bounds(1,:) ) * &
-                  ( [ i, j, k ] + grid_offset ) / real(grid,real32)
-             smallest_bond = modu(get_min_dist(&
-                  basis, tmpvector, .false., &
-                  ignore_list = atom_ignore_list))
-             if( smallest_bond .gt. best_location_bond ) then
-                best_location_bond = smallest_bond
-                best_location = tmpvector
-             end if
-          end do
-       end do
+    do i = 1, size(points,2)
+       smallest_bond = modu(get_min_dist(&
+            basis, [ points(1:3,i) ], .false., &
+            ignore_list = atom_ignore_list))
+       if( smallest_bond .gt. best_location_bond ) then
+          best_location_bond = smallest_bond
+          best_location = points(1:3,i)
+       end if
     end do
 
 
diff --git a/src/fortran/lib/mod_viability.f90 b/src/fortran/lib/mod_viability.f90
index b3526123..c4685d5d 100644
--- a/src/fortran/lib/mod_viability.f90
+++ b/src/fortran/lib/mod_viability.f90
@@ -78,8 +78,8 @@ function get_gridpoints_and_viability(distribs_container, grid, bounds, &
                   ( [ i, j, k ] + grid_offset ) / real(grid,real32)
              do is = 1, basis%nspec
                 atom_loop: do ia = 1, basis%spec(is)%num
-                   do l = 1, size(atom_ignore_list,dim=1), 1
-                      if(all(atom_ignore_list(l,:).eq.[is,ia])) cycle atom_loop
+                   do l = 1, size(atom_ignore_list,dim=2), 1
+                      if(all(atom_ignore_list(:,l).eq.[is,ia])) cycle atom_loop
                    end do
                    if( &
                         get_min_dist_between_point_and_atom( &
diff --git a/src/raffle/raffle.py b/src/raffle/raffle.py
index 1325a881..bb0f57a7 100644
--- a/src/raffle/raffle.py
+++ b/src/raffle/raffle.py
@@ -1503,7 +1503,8 @@ def generate(
                 self, num_structures, stoichiometry,
                 method_ratio={"void": 0.0, "rand": 0.0, "walk": 0.0, "grow": 0.0, "min": 0.0},
                 method_probab=None,
-                seed=None, settings_out_file=None, verbose=0
+                seed=None, settings_out_file=None, verbose=0,
+                calc=None
         ):
             """
             Generate structures using the RAFFLE method.
@@ -1519,12 +1520,23 @@ def generate(
                     DEPRECATED - The ratio of using each method to generate a structure.
                 seed (int):
                     The seed for the random number generator.
-                print_settings (bool):
-                    Boolean whether to print the settings for the generation.
+                settings_out_file (str):
+                    The file to write the settings to.
                 verbose (int):
                     The verbosity level for the generation.
+                calc (ASE calculator):
+                    The calculator to use for the generated structures.
+                
+            Returns:
+                structures (geom_rw.basis_array):
+                    The generated structures.
+                exit_code (int):
+                    The exit code from the generation.
             """
 
+            exit_code = 0
+            structures = None
+
             # check if method_ratio is provided, if so, use it only if method_ratio is not provided
             if method_probab is not None and method_ratio == {"void": 0.0, "rand": 0.0, "walk": 0.0, "grow": 0.0, "min": 0.0}:
                 method_ratio = method_probab
@@ -1548,22 +1560,16 @@ def generate(
             if isinstance(stoichiometry, dict):
                 stoichiometry = Generator.stoichiometry_array(dict=stoichiometry)
 
-            if seed is not None:
-                _raffle.f90wrap_generator__generate__binding__rgt(
-                    this=self._handle,
-                    num_structures=num_structures,
-                    stoichiometry=stoichiometry._handle,
-                    method_ratio=method_ratio_list,
-                    settings_out_file=settings_out_file,
-                    seed=seed, verbose=verbose)
-            else:
-                _raffle.f90wrap_generator__generate__binding__rgt(
-                    this=self._handle,
-                    num_structures=num_structures,
-                    stoichiometry=stoichiometry._handle,
-                    method_ratio=method_ratio_list,
-                    settings_out_file=settings_out_file,
-                    verbose=verbose)
+            exit_code = _raffle.f90wrap_generator__generate__binding__rgt(
+                this=self._handle,
+                num_structures=num_structures,
+                stoichiometry=stoichiometry._handle,
+                method_ratio=method_ratio_list,
+                settings_out_file=settings_out_file,
+                seed=seed, verbose=verbose)
+                
+            structures = self.get_structures(self, calc)[-num_structures:]
+            return structures, exit_code
 
         def get_structures(self, calculator=None):
             """
diff --git a/src/wrapper/f90wrap_mod_generator.f90 b/src/wrapper/f90wrap_mod_generator.f90
index 9c802780..1a77ef17 100644
--- a/src/wrapper/f90wrap_mod_generator.f90
+++ b/src/wrapper/f90wrap_mod_generator.f90
@@ -821,7 +821,9 @@ end subroutine f90wrap_generator__reset_bounds__binding__rgt
 
 subroutine f90wrap_generator__generate__binding__rgt( &
        this, num_structures, stoichiometry, &
-    method_ratio, seed, settings_out_file, verbose)
+       method_ratio, seed, settings_out_file, verbose, exit_code &
+)
+    use raffle__geom_rw, only: basis_type
     use raffle__generator, only: raffle_generator_type, stoichiometry_type
     implicit none
     
@@ -846,6 +848,7 @@ subroutine f90wrap_generator__generate__binding__rgt( &
     character*(*), intent(in), optional :: settings_out_file
     integer, intent(in), optional :: seed
     integer, intent(in), optional :: verbose
+    integer, intent(out), optional :: exit_code
 
     this_ptr = transfer(this, this_ptr)
     stoichiometry_ptr = transfer(stoichiometry, stoichiometry_ptr)
@@ -855,7 +858,8 @@ subroutine f90wrap_generator__generate__binding__rgt( &
          method_ratio=method_ratio, &
          seed=seed, &
          settings_out_file=settings_out_file, &
-         verbose=verbose &
+         verbose=verbose, &
+         exit_code=exit_code &
     )
 end subroutine f90wrap_generator__generate__binding__rgt
 
diff --git a/test/test_evaluator_BTO.f90 b/test/test_evaluator_BTO.f90
index 45d114cb..f5988409 100644
--- a/test/test_evaluator_BTO.f90
+++ b/test/test_evaluator_BTO.f90
@@ -164,12 +164,12 @@ program test_evaluator_BTO
   basis_host%lat(3,:) = [0.0, 0.0, 8.00]
 
   ! set up atom_ignore_list
-  allocate(atom_ignore_list(5,2))
-  atom_ignore_list(1,:) = [1, 2]
-  atom_ignore_list(2,:) = [3, 4]
-  atom_ignore_list(3,:) = [3, 5]
-  atom_ignore_list(4,:) = [3, 6]
-  atom_ignore_list(5,:) = [2, 2]
+  allocate(atom_ignore_list(2,5))
+  atom_ignore_list(:,1) = [1, 2]
+  atom_ignore_list(:,2) = [3, 4]
+  atom_ignore_list(:,3) = [3, 5]
+  atom_ignore_list(:,4) = [3, 6]
+  atom_ignore_list(:,5) = [2, 2]
 
   call basis_host%create_images( &
        max_bondlength = max_bondlength, &
@@ -236,8 +236,8 @@ program test_evaluator_BTO
      write(unit,fmt) basis_host%spec(:)%name
      do is = 1, basis_host%nspec
         atom_loop: do ia = 1, basis_host%spec(is)%num
-           do i = 1, size(atom_ignore_list,1)
-              if( all(atom_ignore_list(i,:).eq.[is,ia]) ) cycle atom_loop
+           do i = 1, size(atom_ignore_list,2)
+              if( all(atom_ignore_list(:,i).eq.[is,ia]) ) cycle atom_loop
            end do
            write(unit,*) basis_host%spec(is)%atom(ia,:3)
         end do atom_loop
@@ -257,29 +257,29 @@ program test_evaluator_BTO
   ! call evaluator
   !-----------------------------------------------------------------------------
   allocate(viability_grid(basis_host%nspec,size(gridpoints,2)))
-  do ia = 1, size(atom_ignore_list,1)
+  do ia = 1, size(atom_ignore_list,2)
      viability_grid(:,:) = 0._real32
      do is = 1, basis_host%nspec
         do i = 1, size(gridpoints,dim=2)
            viability_grid(is,i) = evaluate_point( generator%distributions, &
                 gridpoints(1:3,i), is, &
                 basis_host, &
-                atom_ignore_list(ia:,:), &
+                atom_ignore_list(:,ia:), &
                 [ generator%distributions%bond_info(:)%radius_covalent ] &
            )
         end do
      end do
-     best_loc = maxloc(viability_grid(atom_ignore_list(ia,1),:),dim=1)
+     best_loc = maxloc(viability_grid(atom_ignore_list(1,ia),:),dim=1)
      ltmp1 = .false.
-     if(any(viability_grid(atom_ignore_list(ia,1),:) .gt. 0._real32) )then
-        do ja = ia, size(atom_ignore_list,1), 1
-           if( atom_ignore_list(ja,1) .ne. atom_ignore_list(ia,1) ) cycle
+     if(any(viability_grid(atom_ignore_list(1,ia),:) .gt. 0._real32) )then
+        do ja = ia, size(atom_ignore_list,2), 1
+           if( atom_ignore_list(1,ja) .ne. atom_ignore_list(1,ia) ) cycle
            if( &
                 all( &
                      abs( &
                           gridpoints(1:3,best_loc) - &
-                          basis_host%spec(atom_ignore_list(ja,1))%atom( &
-                               atom_ignore_list(ja,2),:3 &
+                          basis_host%spec(atom_ignore_list(1,ja))%atom( &
+                               atom_ignore_list(2,ja),:3 &
                           ) &
                      ) .lt. tolerance + 1.E-6_real32 &
                 ) &
@@ -292,15 +292,15 @@ program test_evaluator_BTO
           success &
      )
      if( .not. ltmp1 ) then
-        write(*,*) "species: ", atom_ignore_list(ia,1)
+        write(*,*) "species: ", atom_ignore_list(1,ia)
         write(*,*) "Best location: ", best_loc
         write(*,*) "viability: ", &
-             viability_grid(atom_ignore_list(ia,1),best_loc)
+             viability_grid(atom_ignore_list(1,ia),best_loc)
         write(*,*) "Gridpoint: ", gridpoints(1:3,best_loc)
      end if
      call basis_host%update_images( &
           max_bondlength = max_bondlength, &
-          is = atom_ignore_list(ia,1), ia = atom_ignore_list(ia,2) &
+          is = atom_ignore_list(1,ia), ia = atom_ignore_list(2,ia) &
      )
   end do
   close(unit)
diff --git a/test/test_evaluator_C.f90 b/test/test_evaluator_C.f90
index 9cf51b9e..5b510ae3 100644
--- a/test/test_evaluator_C.f90
+++ b/test/test_evaluator_C.f90
@@ -149,11 +149,11 @@ program test_evaluator
   basis_host%lat(3,:) = [0.0, 0.0, 7.121490218]
 
   ! set up atom_ignore_list
-  allocate(atom_ignore_list(8,2))
-  do i = 1, size(atom_ignore_list,1)
-     atom_ignore_list(i,1) = 1
-     atom_ignore_list(i,2) = &
-          basis_host%spec(1)%num - size(atom_ignore_list,1) + i
+  allocate(atom_ignore_list(2,8))
+  do i = 1, size(atom_ignore_list,2)
+     atom_ignore_list(1,i) = 1
+     atom_ignore_list(2,i) = &
+          basis_host%spec(1)%num - size(atom_ignore_list,2) + i
   end do
 
   call basis_host%create_images( &
@@ -220,8 +220,8 @@ program test_evaluator
      write(unit,fmt) basis_host%spec(:)%name
      do is = 1, basis_host%nspec
         atom_loop: do ia = 1, basis_host%spec(is)%num
-           do i = 1, size(atom_ignore_list,1)
-              if( all(atom_ignore_list(i,:).eq.[is,ia]) ) cycle atom_loop
+           do i = 1, size(atom_ignore_list,2)
+              if( all(atom_ignore_list(:,i).eq.[is,ia]) ) cycle atom_loop
            end do
            write(unit,*) basis_host%spec(is)%atom(ia,:3)
         end do atom_loop
@@ -241,24 +241,24 @@ program test_evaluator
   ! call evaluator
   !-----------------------------------------------------------------------------
   allocate(viability_grid(basis_host%nspec,size(gridpoints,2)))
-  do ia = 1, size(atom_ignore_list,1)
+  do ia = 1, size(atom_ignore_list,2)
      viability_grid(:,:) = 0._real32
      do i = 1, size(gridpoints,dim=2)
         viability_grid(1,i) = evaluate_point( generator%distributions, &
-             gridpoints(1:3,i), atom_ignore_list(ia,1), basis_host, &
-             atom_ignore_list(ia:,:), &
+             gridpoints(1:3,i), atom_ignore_list(1,ia), basis_host, &
+             atom_ignore_list(:,ia:), &
              [ generator%distributions%bond_info(:)%radius_covalent ] &
         )
      end do
-     best_loc = maxloc(viability_grid(atom_ignore_list(ia,1),:),dim=1)
+     best_loc = maxloc(viability_grid(atom_ignore_list(1,ia),:),dim=1)
      ! Check point is correct
      ltmp1 = .false.
-     do ja = ia, size(atom_ignore_list,1), 1
+     do ja = ia, size(atom_ignore_list,2), 1
         if( &
              all( &
                   abs( &
                        gridpoints(1:3,best_loc) - &
-                       basis_host%spec(1)%atom(atom_ignore_list(ja,2),:3) &
+                       basis_host%spec(1)%atom(atom_ignore_list(2,ja),:3) &
                   ) .lt. tolerance + 1.E-6_real32 &
              ) &
         ) ltmp1 = .true.
@@ -270,7 +270,7 @@ program test_evaluator
      )
      call basis_host%update_images( &
           max_bondlength = max_bondlength, &
-          is = 1, ia = atom_ignore_list(ia,2) &
+          is = 1, ia = atom_ignore_list(2,ia) &
      )
   end do
 
diff --git a/test/test_geom_extd.f90 b/test/test_geom_extd.f90
index 831b50fc..f3e4319f 100644
--- a/test/test_geom_extd.f90
+++ b/test/test_geom_extd.f90
@@ -88,10 +88,10 @@ subroutine test_update_iamges(basis, success)
     logical, intent(inout) :: success
 
     type(extended_basis_type) :: basis_copy
-    integer, dimension(1,2) :: atom_ignore_list
+    integer, dimension(2,1) :: atom_ignore_list
 
     call basis_copy%copy(basis)
-    atom_ignore_list(1,:) = [1, 1]
+    atom_ignore_list(:,1) = [1, 1]
 
     ! Create images
     call basis_copy%create_images( &
diff --git a/test/test_place_methods.f90 b/test/test_place_methods.f90
index f5d292f4..10987526 100644
--- a/test/test_place_methods.f90
+++ b/test/test_place_methods.f90
@@ -67,12 +67,6 @@ program test_place_methods
   basis%lat(3,3) = 5.0_real32
 
 
-  !-----------------------------------------------------------------------------
-  ! test place_method_void
-  !-----------------------------------------------------------------------------
-  call test_place_method_void(basis, success)
-
-
   !-----------------------------------------------------------------------------
   ! set up distribution functions
   !-----------------------------------------------------------------------------
@@ -95,8 +89,8 @@ program test_place_methods
   call generator%distributions%host_system%set_element_map( &
        generator%distributions%element_info &
   )
-  allocate(atom_ignore_list(1, 2))
-  atom_ignore_list(1,:) = [1,2]
+  allocate(atom_ignore_list(2,1))
+  atom_ignore_list(:,1) = [1,2]
   seed = 0
   call random_seed(size=num_seed)
   allocate(seed_arr(num_seed))
@@ -109,6 +103,12 @@ program test_place_methods
   )
 
 
+  !-----------------------------------------------------------------------------
+  ! test place_method_void
+  !-----------------------------------------------------------------------------
+  call test_place_method_void(generator%distributions, basis, success)
+
+
   !-----------------------------------------------------------------------------
   ! test place_method_rand
   !-----------------------------------------------------------------------------
@@ -199,10 +199,12 @@ program test_place_methods
 
 contains
 
-  subroutine test_place_method_void(basis, success)
+  subroutine test_place_method_void(distributions, basis, success)
+    use raffle__viability, only: get_gridpoints_and_viability
     implicit none
     logical, intent(inout) :: success
     type(basis_type), intent(in) :: basis
+    type(distribs_container_type), intent(in) :: distributions
 
     integer :: i
     type(extended_basis_type) :: basis_copy
@@ -213,21 +215,40 @@ subroutine test_place_method_void(basis, success)
     integer, dimension(:,:), allocatable :: atom_ignore_list
     real(real32), dimension(3) :: tolerance
     real(real32), dimension(2,3) :: bounds
+    real(real32), dimension(:,:), allocatable :: gridpoints, viability_grid
 
     ! Initialise test data
     grid = [10, 10, 10]
     bounds(1,:) = 0.0_real32
     bounds(2,:) = 1.0_real32
-    allocate(atom_ignore_list(1, 2))  ! No atoms to ignore
-    atom_ignore_list(1,:) = [1,2]
+    allocate(atom_ignore_list(2,1))  ! No atoms to ignore
+    atom_ignore_list(:,1) = [1,2]
     grid_offset = [0.5_real32, 0.5_real32, 0.5_real32]
 
     ! Initialise basis
     call basis_copy%copy(basis)
+    call basis_copy%create_images( &
+         max_bondlength = 6._real32, &
+         atom_ignore_list = atom_ignore_list &
+    )
+
+    !---------------------------------------------------------------------------
+    ! set up gridpoints
+    !---------------------------------------------------------------------------
+    gridpoints = get_gridpoints_and_viability( &
+         distributions, &
+         grid, &
+         bounds, &
+         basis_copy, &
+         [ 1 ], &
+         [ distributions%bond_info(:)%radius_covalent ], &
+         atom_ignore_list, &
+         grid_offset = grid_offset &
+    )
 
     ! Call the void subroutine
     point = place_method_void( &
-         grid, grid_offset, bounds, basis_copy, &
+         gridpoints, basis_copy, &
          atom_ignore_list, &
          viable &
     )
diff --git a/test/test_viability.f90 b/test/test_viability.f90
index f340ea9c..f49210d4 100644
--- a/test/test_viability.f90
+++ b/test/test_viability.f90
@@ -63,8 +63,8 @@ subroutine test_get_gridpoints_and_viability(basis, success)
     grid = [10, 10, 10]
     bounds(1,:) = 0.0_real32
     bounds(2,:) = 1.0_real32
-    allocate(atom_ignore_list(1, 2))  ! No atoms to ignore
-    atom_ignore_list(1,:) = [1,2]
+    allocate(atom_ignore_list(2,1))  ! No atoms to ignore
+    atom_ignore_list(:,1) = [1,2]
     allocate(radius_list(1))
     radius_list = 1.0_real32
     lowtol = 0.5_real32
@@ -133,8 +133,8 @@ subroutine test_update_gridpoints_and_viability(basis, success)
     grid = [10, 10, 10]
     bounds(1,:) = 0.0_real32
     bounds(2,:) = 1.0_real32
-    allocate(atom_ignore_list(1, 2))  ! No atoms to ignore
-    atom_ignore_list(1,:) = [1,2]
+    allocate(atom_ignore_list(2,1))  ! No atoms to ignore
+    atom_ignore_list(:,1) = [1,2]
     allocate(radius_list(1))
     radius_list = 1.0_real32 !!! NO!!! USING CARBON RADIUS
     lowtol = 0.5_real32