diff --git a/.gitignore b/.gitignore
index 06de34a..d0f861d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,9 @@
 __pycache__
 inout_files/
+inout_files_control/
+inout_files_heli_only/
+inout_files_heli_shield/
+inout_files_shield_only/
 *.txt
 *.h5
 *.log
diff --git a/README.md b/README.md
index 42d0be4..4d3982b 100644
--- a/README.md
+++ b/README.md
@@ -1,38 +1,77 @@
 # HeliMax
-This repository aims to show how to generate gprMax input files using NumPy arrays, and how to perform simulation using these input files. The models generated and performed are inspired by an airborne ice Radar setup from ETH Zurich. This work takes place is a master thesis project. The ultimate goal of these simulations is to determine if a shielding between the GPR and the helicopter would be effective and if so, what geometry and properties it should have. The framework can be used to perform any type of layer-based models.
- 
-Public description of the project : https://vaw.ethz.ch/en/research/glaciology/research-projects.html
- 
-# The setup
-The setup is composed of multiple python files, which the most important one is `main.py`. The folder `inout_files` stores  input and output files, as well as the figures. This repository should be downloaded and placed into the gprMax folder. For example :
- 
-`home/user/gprMax/HeliMax`
- 
- You must have a gprMax environment correctly installed : https://github.com/gprMax/gprMax
- You must install additionnal requirments (´conda requirment.txt´)
- 
-# How it works
-The script is executed using the terminal. `python main.py` generates the model, you can check your geometry and the files generated. `python main.py --run` regenerate all the files and launch the gprMax simulation.
-
-## Generate models
-The models can be generated using simple NumPy arrays. Any kind of geometry can be generated, on few conditions:
-1. The array is composed of integer (not physical properties)
-2. Every integer correspond to a specific material
-3. All materials are described in a `.txt` file (see the file for the structure)
-4. The order of the materials matters
-5. The matrix of integer is stored in a `.h5` file
-6. The `.in` file contains all the information, plus the `.txt` and `.h5` file paths
-
-## gprMax simulations
-The documentation of gprMax give some simple commands to run simulations. For example:
- 
-`python -m gprMax user_models/cylinder_Ascan_2D.in`, 
- 
-where `cylinder_Ascan.in` is a basic input file from gprMax. In certain cases, it's an advantage to use a GPU to run gprMax simulations. This can be done by using the `-gpu` argument.
-
-# Multi-GPU clusters
-For multi-GPU users, the installation requires some attention. In the same conda environment used for gprMax, some packages must be installed. The three main packages to install are `pycuda`, `openmpy` and `mpi4py`. This must be done using:
-
-`conda install -c conda-forge pycuda mpi4py openmpi`
-
-Then, in certain cases, it will be necessary to remove tha package `mpich` which is installed by conda when installing the package `mpi4py`. `pip` is not used as it will struggle to build wheels to install the packages. 
+
+This repository contains a Python program designed to run gprMax simulations specifically tailored for glaciology research. It utilizes ground-penetrating radar (GPR) simulation to understand and analyze sub-glacial features, including bedrock roughness and ice thickness. The program supports model creation, simulation execution, file merging, and result plotting, with the capability to specify GPU usage for computational efficiency.
+
+## Features
+
+- **Model Generation**: Dynamically generates simulation models based on user-defined parameters.
+- **Simulation Execution**: Runs simulations using specified GPUs to leverage computational resources effectively.
+- **File Merging**: Merges simulation output files for consolidated analysis.
+- **Result Plotting**: Plots the simulation results for visual analysis.
+
+## Requirements
+
+- Python 3.x
+- gprMax
+- tqdm
+- multiprocessing
+
+Please ensure you have gprMax installed and configured correctly on your system. You can find installation instructions for gprMax at [gprMax's GitHub page](https://github.com/gprMax/gprMax).
+
+## Installation
+
+Clone this repository to your local machine, make sure that it's in you folder `gprMax`:
+
+```bash
+git clone https://github.com/bart-inho/helimax.git
+cd helimax
+```
+
+Ensure you have the necessary Python dependencies installed:
+
+```bash
+pip install -r requirements.txt
+```
+
+## Usage
+
+The program is executed via the command line, offering several options to control its operation:
+
+- `--run`: Execute the simulation.
+- `--model`: Generate the simulation model.
+- `--merge`: Merge simulation output files.
+- `--plot`: Plot the simulation results.
+- `--name`: Specify the path for input and output files. This argument is required.
+- `--gpus`: Specify comma-separated GPU IDs to use for the simulation (e.g., `0,1,2,3`).
+
+### Examples
+
+Generate a model:
+
+```bash
+python main.py --model --name path/to/files
+```
+
+Run simulations using GPUs 0 and 1:
+
+```bash
+python main.py --run --name path/to/files --gpus 0,1
+```
+
+Merge simulation files and plot results:
+
+```bash
+python main.py --merge --plot --name path/to/files
+```
+
+## Customization
+
+You can customize the simulation parameters such as antenna spacing, measurement number, and dimensions directly within the script. For more advanced configurations, edit the model generation logic as needed.
+
+## Contributing
+
+Contributions to improve the program or address issues are welcome. Please feel free to submit a pull request or open an issue.
+
+## Acknowledgments
+
+This project utilizes [gprMax](https://github.com/gprMax/gprMax), an open-source project for simulating GPR. We acknowledge the developers and contributors of gprMax for their valuable work.
\ No newline at end of file
diff --git a/conda_requirment.txt b/conda_requirment.txt
deleted file mode 100644
index 579f9f3..0000000
--- a/conda_requirment.txt
+++ /dev/null
@@ -1,7 +0,0 @@
-# Install dependencies
-mpi4py
-openmpi
-cuda-nvcc
-
-# Set environment variables
-export CPLUS_INCLUDE_PATH=/usr/local/cuda/include:${CPLUS_INCLUDE_PATH}
\ No newline at end of file
diff --git a/main.py b/main.py
index 5f3d9c9..c23817e 100644
--- a/main.py
+++ b/main.py
@@ -1,91 +1,87 @@
 # This program is used to run the gprMax simulation applied to glaciology.
-from services.file_service import FileService
 from services.folder_init import InitializeFolders
-from models.simulation_model import SimulationModel
-from models.material import Material
+from models.model_generation_logic import ModelGenerationLogic
 from simulations.simulation_runner import SimulationRunner
-from simulations.simulation_plot_profile import PlotProfile  
+from simulations.simulation_plot_profile import PlotProfile
 import argparse
+from tqdm import tqdm  # Import tqdm
+from concurrent.futures import ThreadPoolExecutor
+from multiprocessing import Pool
 
 def main():
     parser = argparse.ArgumentParser()
     parser.add_argument('--run', action='store_true', help='run the simulation')
+    parser.add_argument('--model', action='store_true', help='create the model')
+    parser.add_argument('--merge', action='store_true', help='merge the simulation files')
     parser.add_argument('--plot', action='store_true', help='plot the simulation')
     parser.add_argument('--rough', action='store_true', help='rough bedrock')
+    parser.add_argument('--name', type=str, help='Path to input and output files', required=True)
+    parser.add_argument('--gpus', type=str, help='Comma-separated list of GPU IDs to use', default='0')
+
     args = parser.parse_args()
 
     # Initialize folders
-    InitializeFolders.check_and_create_directories()
-
-    # Initialize Materials
-    # Change the material names in the "Material" class 
-    freespace = Material(1. , 0.   , 1., 0., 'freespace') # Free space
-    glacier   = Material(3.2, 5.e-8, 1., 0., 'glacier'  ) # Glacier
-    bedrock   = Material(5. , 1.e-2, 1., 0., 'bedrock'  ) # Bedrock
-    metal     = Material(1. , 'inf', 1., 0., 'metal'    ) # Helico
+    model_name    = 'model'
+    path_to_files = args.name
+    InitializeFolders.check_and_create_directories(path_to_files)
     
-    # Initialize SimulationModel
-    model_name    = 'test_roughness'
-    inout_files   = 'inout_files/'
-    path_to_files = inout_files + model_name
+    antenna_spacing = 4  # Change antenna spacing in [m] here
 
-    # Generate model
-    model = SimulationModel(model_name, 
-                            45, 10, 100, 
-                            [0.08, 0.08, 0.08], # Change discretisation if needed here
-                            [freespace, glacier, bedrock, metal], # Change name of materials here
-                            inout_files)
-    
-    # Generate base model
-    model.generate_base_glacier()
-    model.generate_curved_bedrock_glacier([-10, 5, -200], # center of the curvature [m]
-                             100,            # radius of the curvature [m]
-                             args.rough)
+    dis = 0.06
+    measurement_number = 200 # number of traces
 
-    measurement_number = 24 # number of traces
-    antenna_spacing    = 4  # Change antenna spacing in [m] here
-    measurement_step   = model.calculate_measurment_step(measurement_number, 
-                                                         antenna_spacing) # Change antenna spacing in m here
-     
-    # Add antenna positions
-    transceiver1 = [round(25 * model.discrete[0]), # 25 cells of buffer (20 minimum)    
-                    round(model.y_size/2        ),
-                    round(model.z_size/5-.5     )]
-    
-    receiver1    = [round(25 * model.discrete[0] + antenna_spacing), # 25 cells of buffer (20 minimum)
-                    round(model.y_size/2                          ),
-                    round(model.z_size/5-.5                       )]
-        
-    #Plot initial model
-    model.plot_initial_model(transceiver1, receiver1)
+    x_m = 60
+    y_m = 10
+    z_m = 100
 
-    # Call FileService to write files
-    FileService.write_materials_file(model.path + model.name + '_materials', 
-                                     model.materials)
-    
-    FileService.write_h5_file(model.path + model.name + '_h5', 
-                              model)
+    measurement_step = (x_m-7.5) / measurement_number
+    print(measurement_step)
 
-    FileService.write_input_file(model, 
-                                path_to_files, 
-                                path_to_files + '_materials', 
-                                path_to_files + '_h5', 
-                                25e6,   # Change frequency in Hz here
-                                transceiver1, receiver1, 
-                                measurement_step, 
-                                1000e-9) # Change time window in s here
-        
-    # Run simulation
+    antenna_x = round(30 * dis) # 30 cells of buffer (20 minimum
+    antenna_y = round(y_m/2)
+    antenna_z = round(20)
+
+    rope_length = 15 # length of the rope in [m]
+
+    if args.model:
+        with ThreadPoolExecutor(max_workers=18) as executor:
+            import concurrent.futures
+            futures = [executor.submit(ModelGenerationLogic.model_generation_logic, idx, 
+                                        model_name, 
+                                        x_m, y_m, z_m, 
+                                        dis, antenna_spacing, 
+                                        rope_length, 
+                                        path_to_files, 
+                                        antenna_x+(idx*measurement_step), antenna_y, antenna_z,
+                                        measurement_step,
+                                        args) for idx in range(measurement_number)]
+            
+            # Process the futures as they complete (to maintain order, if necessary)
+            for future in tqdm(concurrent.futures.as_completed(futures), total=measurement_number, desc='Creating Models'):
+                result = future.result()  # result returned from create_model
+                
     if args.run:
-        simulation_runner = SimulationRunner(model)
-        simulation_runner.run_simulation(measurement_number)
-        simulation_runner.merge_files(True)
-        
-    # Plot profile
+        gpu_ids = [int(gpu) for gpu in args.gpus.split(',')]
+
+        # Prepare arguments for each simulation
+        total_gpus = len(gpu_ids)
+        simulation_args = [(args.name, 'model', idx, gpu_ids[idx % total_gpus]) for idx in range(200)]
+
+        # Use multiprocessing Pool to run simulations in parallel
+        with Pool(processes=total_gpus) as pool:
+            list(tqdm(pool.imap(SimulationRunner.run_simulation, simulation_args), total=len(simulation_args)))
+
+        print("All simulations completed.")
+
+    if args.merge:
+        # Merge the files
+        SimulationRunner.merge_files(False, path_to_files, model_name)
+
     if args.plot:
-        plot_profile = PlotProfile(model.path + model.name + '_merged.out', 'Ey')
+        # Plot the simulation
+        plot_profile = PlotProfile(path_to_files + model_name + '_merged.out', 'Ey')
         plot_profile.get_output_data()
         plot_profile.plot()
-
+        
 if __name__ == "__main__":
     main()
\ No newline at end of file
diff --git a/models/model_generation_logic.py b/models/model_generation_logic.py
new file mode 100644
index 0000000..166651e
--- /dev/null
+++ b/models/model_generation_logic.py
@@ -0,0 +1,83 @@
+from services.file_service import FileService
+from models.simulation_model import SimulationModel
+from models.material import Material
+
+class ModelGenerationLogic:
+    def model_generation_logic(idx, 
+                    model_name, 
+                    x_m, y_m, z_m, dis, 
+                    antenna_spacing, 
+                    rope_length, 
+                    path_to_files, 
+                    antenna_x, antenna_y, antenna_z,
+                    measurement_step,
+                    args):
+
+        # Initialize Materials
+        # Change the material names in the "Material" class 
+        freespace = Material(1. , 0.   , 1., 0., 'freespace') # Free space
+        glacier   = Material(3.2, 5.e-8, 1., 0., 'glacier'  ) # Glacier
+        bedrock   = Material(5. , 1.e-2, 1., 0., 'bedrock'  ) # Bedrock
+        metal     = Material(1. , 'inf', 1., 0., 'metal'    ) # Helico
+        shield    = Material(6. , 1.e-2, 1., 0, 'shield'    ) # Shield
+
+        # Initialize SimulationModel
+        if idx == 0:
+            model_name_it = model_name
+        elif idx > 0:
+            model_name_it = model_name + str(idx)
+
+        # Generate model
+        model = SimulationModel(model_name_it, 
+                                x_m, y_m, z_m, 
+                                [dis, dis, dis], # Change discretisation if needed here
+                                [freespace, glacier, bedrock, metal, shield], # Change name of materials here
+                                path_to_files)
+        
+        # Generate base model
+        model.generate_base_glacier()
+        model.generate_curved_bedrock_glacier([-10, 5, -200], # center of the curvature [m]
+                                            100,            # radius of the curvature [m]
+                                            args.rough)
+        
+        model.add_CI_shield(antenna_x, antenna_y, antenna_z, dis, antenna_spacing, dis)
+        # model.add_DMC_Tx(antenna_x, antenna_y, antenna_z, dis, antenna_spacing, dis)
+        # model.add_DMC_Rx(antenna_x, antenna_y, antenna_z, dis, antenna_spacing, dis)
+
+        # model.add_3D_oval_shape([antenna_x , antenna_y, antenna_z - rope_length], # center of the curvature [m]
+        #                         [10, .5, .5],            # radius of the curvature [m]])
+        #                         dis*20)
+        # model.add_3D_oval_shape([antenna_x + 3.5, antenna_y, antenna_z - rope_length - 1.5], # center of the curvature [m]
+        #                         [20, .1, .1],            # radius of the curvature [m]])
+        #                         dis*100)
+        # model.add_3D_oval_shape([antenna_x + antenna_spacing - 1, antenna_y, antenna_z - rope_length], # center of the curvature [m]
+        #                         [5, 2, 2],            # radius of the curvature [m]])
+        #                         dis*3)
+        
+        transceiver1 = [antenna_x, # 25 cells of buffer (20 minimum)    
+                        antenna_y,
+                        antenna_z]
+        
+        receiver1    = [antenna_x + antenna_spacing, # 25 cells of buffer (20 minimum)
+                        antenna_y,
+                        antenna_z]
+            
+        #Plot initial model each 20 steps
+        # if idx%20 == 0:
+            # model.plot_initial_model(transceiver1, receiver1)
+
+        # Call FileService to write files
+        FileService.write_materials_file(model.path + model.name + '_materials', 
+                                        model.materials)
+        
+        FileService.write_h5_file(model.path + model.name + '_h5', 
+                                model)
+
+        FileService.write_input_file(model, 
+                                    path_to_files+model_name_it, 
+                                    path_to_files+model_name_it + '_materials', 
+                                    path_to_files+model_name_it + '_h5', 
+                                    25e6,   # Change frequency in Hz here
+                                    transceiver1, receiver1, 
+                                    measurement_step, 
+                                    1000e-9) # Change time window in s here
\ No newline at end of file
diff --git a/models/simulation_model.py b/models/simulation_model.py
index f41095e..7a36974 100644
--- a/models/simulation_model.py
+++ b/models/simulation_model.py
@@ -1,5 +1,6 @@
 import numpy as np
 import matplotlib.pyplot as plt
+import matplotlib.colors as mcolors
 
 class SimulationModel:
     # This class is used to store the model information and generate the base model.
@@ -23,10 +24,69 @@ def generate_base_glacier(self):
         nx = int(self.x_size / self.discrete[0])
         ny = int(self.y_size / self.discrete[1])
         nz = int(self.z_size / self.discrete[2])
-        
+
         self.model = np.zeros((nx, ny, nz)) # Free space = 0
-        self.model[:, :, round(nz/3):nz] = 1 # Glacier = 1
-        # self.model[:, :, 0:10] = 3 # Metal = 3
+        self.model[:, :, round(35/self.discrete[2]):nz] = 1 # Glacier = 1
+
+    def add_CI_shield(self, antenna_x, antenna_y, antenna_z, thick, antenna_spacing, dis):
+        
+        right_corner = round((antenna_x)/dis)
+        left_corner = round((antenna_x + antenna_spacing + .1)/dis)
+        # right_corner = round((antenna_x - 2)/dis)
+        # left_corner = round((antenna_x + antenna_spacing + 2)/dis)
+        height = round((antenna_z-2)/dis)
+        thickness = round(thick/dis)
+        front = round((antenna_y - 2)/dis)
+        back = round((antenna_y + 2)/dis)
+        # front = round((antenna_y - 3)/dis)
+        # back = round((antenna_y + 3)/dis)
+        self.model[right_corner:left_corner, front:back, height-thickness:height] = 3
+
+    def add_DMC_Tx(self, antenna_x, antenna_y, antenna_z, thick, antenna_spacing, dis):
+        
+        right_corner = round((antenna_x - 2*dis)/dis)
+        left_corner = round((antenna_x + 4*dis)/dis)
+        height = round((antenna_z-dis)/dis)
+        thickness = round(thick/dis)
+        front = round((antenna_y - 2)/dis)
+        back = round((antenna_y + 2)/dis)
+        self.model[right_corner:left_corner, front:back, height-thickness:height] = 3
+
+    def add_DMC_Rx(self, antenna_x, antenna_y, antenna_z, thick, antenna_spacing, dis):
+        
+        right_corner = round((antenna_x + antenna_spacing - 2*dis)/dis)
+        left_corner = round((antenna_x + antenna_spacing + 4*dis)/dis)
+        height = round((antenna_z-dis)/dis)
+        thickness = round(thick/dis)
+        front = round((antenna_y - 2)/dis)
+        back = round((antenna_y + 2)/dis)
+        self.model[right_corner:left_corner, front:back, height-thickness:height] = 3
+
+    def add_3D_oval_shape(self, center, radii, thickness):
+        nx, ny, nz = self.model.shape
+        center = np.array(center) / self.discrete
+        radii = np.array(radii) / self.discrete / 2
+        thickness = thickness / self.discrete[0]
+
+        # Create index arrays for i, j, k
+        i, j, k = np.mgrid[:nx, :ny, :nz]
+
+        # Calculate the squared distances from the center divided by the radii squared
+        # to create an ellipsoid equation of the form (x^2/a^2 + y^2/b^2 + z^2/c^2 = 1)
+        ellipsoid = (((i - center[0]) / radii[0])**2 +
+                     ((j - center[1]) / radii[1])**2 +
+                     ((k - center[2]) / radii[2])**2  )
+
+        # Define the outer and inner surfaces of the shell
+        outer_surface = ellipsoid <= 1
+        inner_surface = ellipsoid <= (1 - thickness/radii[0])**2
+
+        # Create shell between the inner and outer surfaces
+        shell_mask = outer_surface & ~inner_surface
+
+        # Update the model with the new values where the shell should be
+        self.model[shell_mask] = 3  # Assign a value for the shell material
+        self.model[inner_surface] = 0  # Assign a value for the inner surface
 
     def generate_curved_bedrock_glacier(self, center, r, arg_rough):
         nx, ny, nz = self.model.shape
@@ -56,19 +116,47 @@ def generate_curved_bedrock_glacier(self, center, r, arg_rough):
         self.model[ice_mask] = 1  # Ice = 1
 
     def plot_initial_model(self, transceiver, receiver):
+        """
+        Plot the initial model
+
+        Parameters:
+        self (SimulationModel): the model to plot
+        transceiver (np.array): the transceiver position
+        receiver (np.array): the receiver position
+
+        Returns:
+        None
+        """
 
-        X, Y = np.meshgrid(np.arange(0, self.x_size, self.discrete[0]), 
-                        np.arange(0, self.z_size, self.discrete[2]))
+        nx = self.model.shape[0]
+        nz = self.model.shape[2]
+
+        model = self.model[:, round(self.y_size/self.discrete[1]/2), :].T
+
+        X, Y = np.meshgrid(np.linspace(0, self.x_size + self.discrete[0], nx+1), 
+                        np.linspace(0, self.z_size + self.discrete[2], nz+1))
         
-        plt.pcolormesh(X, Y, self.model[:, round(transceiver[1]*self.discrete[1]), :].T)
-        plt.scatter(transceiver[0], transceiver[2])
-        plt.scatter(receiver[0], receiver[2])
+        # Plot parameters
+        colors = ['white', 'blue', 'gray', 'black']  # Define colors for each category
+
+        # Create a colormap from the list of colors
+        cmap = mcolors.ListedColormap(colors)
+
+        # Create a boundary norm with the categories as boundaries
+        bounds = [0, 1, 2, 3, 4]  # Set boundaries including the right edge for the last category
+        norm = mcolors.BoundaryNorm(bounds, cmap.N)
+
+        plt.pcolormesh(X, Y, model, cmap=cmap, norm=norm, shading='auto')
+        cbar = plt.colorbar(ticks=[0.5, 1.5, 2.5, 3.5], boundaries=bounds, orientation='horizontal', pad=0.11, aspect=50)  # Set tick positions at the center of each category
+        cbar.set_ticklabels(['Free Space', 'Glacier Ice', 'Bedrock', 'Metal'])  # Label each tick with text
+        plt.scatter(transceiver[0], transceiver[2], .5, color='red', label='Transceiver')  # Mark transceiver
+        plt.scatter(receiver[0], receiver[2], .5, color='green', label='Receiver')  # Mark receiver
+        # plt.legend(loc='upper right')
         plt.gca().invert_yaxis()
-        plt.gca().set_aspect('equal')
-        plt.colorbar()
-        plt.clim(0, 3)
+        plt.gca().set_aspect('equal', adjustable='box')
         plt.ylabel('depth [m]')
         plt.xlabel('distance [m]')
-        plt.title(self.name)
-        plt.savefig(self.path+'/figures/'+self.name+'.png')
-        plt.close()        
\ No newline at end of file
+        plt.title('Moon Lava Tube GPR Test')
+        plt.tight_layout()
+        plt.savefig(self.path+'/figures/'+self.name+'.png', format='png')  # Use efficient format
+        plt.close()
\ No newline at end of file
diff --git a/services/folder_init.py b/services/folder_init.py
index c0b3366..a00c667 100644
--- a/services/folder_init.py
+++ b/services/folder_init.py
@@ -3,6 +3,6 @@
 class InitializeFolders:
 
     @staticmethod
-    def check_and_create_directories():
+    def check_and_create_directories(folder_name):
         # Check if inout_files exists and create if not
-        Path('inout_files/figures').mkdir(parents=True, exist_ok=True)
\ No newline at end of file
+        Path(folder_name+'figures').mkdir(parents=True, exist_ok=True)
diff --git a/simulations/simulation_runner.py b/simulations/simulation_runner.py
index 4ad6ca7..6673e35 100644
--- a/simulations/simulation_runner.py
+++ b/simulations/simulation_runner.py
@@ -1,12 +1,27 @@
 from gprMax.gprMax import api
 from tools.outputfiles_merge  import merge_files
+from queue import Queue
+import os
+from multiprocessing import current_process
+
 
 class SimulationRunner:
-    def __init__(self, simulation_model):
-        self.simulation_model = simulation_model
-
-    def run_simulation(self, measurement_number):
-        api(self.simulation_model.path + self.simulation_model.name + '.in', mpi = 6, gpu = [0, 1, 2, 3, 4, 5], n = measurement_number)
-    
-    def merge_files(self, remove_files):
-        merge_files(self.simulation_model.path + self.simulation_model.name, removefiles = remove_files)
\ No newline at end of file
+    def run_simulation(args):
+        path_to_files, model_name, idx, gpu_id = args
+        model_input_file_name = model_name if idx == 0 else f"{model_name}{idx}"
+        
+        print(f"Process {current_process().name} running on GPU {gpu_id}: {os.path.join(path_to_files, model_input_file_name + '.in')}")
+        api(os.path.join(path_to_files, model_input_file_name + ".in"), mpi=False, gpu=[gpu_id], n = 1)
+
+    def create_gpu_queue(gpu_ids):
+        """
+        Create a queue populated with the specified GPU IDs.
+        """
+        gpu_queue = Queue() # Create a queue
+        for gpu_id in gpu_ids: # Iterate over the specified GPU IDs
+            gpu_queue.put(gpu_id)  # Populate the queue with specified GPU IDs
+        return gpu_queue
+
+
+    def merge_files(remove_files, path_to_files, model_name):
+        merge_files(path_to_files+model_name, removefiles = remove_files)
\ No newline at end of file