diff --git a/netpyne/sim/__init__.py b/netpyne/sim/__init__.py
index d43f0f007..3529f0737 100644
--- a/netpyne/sim/__init__.py
+++ b/netpyne/sim/__init__.py
@@ -66,6 +66,7 @@
     _init_stim_randomizer,
     unique,
     checkMemory,
+    close,
 )
 
 # import utils functions to manipulate objects
diff --git a/netpyne/sim/utils.py b/netpyne/sim/utils.py
index e233a4adc..34454f755 100644
--- a/netpyne/sim/utils.py
+++ b/netpyne/sim/utils.py
@@ -13,6 +13,7 @@
 from time import time
 import hashlib
 import array
+import sys
 from numbers import Number
 from collections import OrderedDict
 from neuron import h  # Import NEURON
@@ -914,7 +915,8 @@ def clearAll():
 
     # clean up
     sim.pc.barrier()
-    sim.pc.gid_clear()  # clear previous gid settings
+    #TODO are the references within containers cleared before the containers are?
+
 
     # clean cells and simData in all nodes
     if hasattr(sim, 'net'):
@@ -963,54 +965,61 @@ def clearAll():
             matplotlib.pyplot.clf()
             matplotlib.pyplot.close('all')
 
-    # clean rxd components
-    if hasattr(sim.net, 'rxd'):
-
-        sim.clearObj(sim.net.rxd)
-
-        if 'rxd' not in globals():
-            try:
-                from neuron import crxd as rxd
-            except:
-                pass
-        # try:
-        for r in rxd.rxd._all_reactions[:]:
-            if r():
-                rxd.rxd._unregister_reaction(r)
-
-        for s in rxd.species._all_species:
-            if s():
-                s().__del__()
-
-        rxd.region._all_regions = []
-        rxd.region._region_count = 0
-        rxd.region._c_region_lookup = None
-        rxd.species._species_counts = 0
-        rxd.section1d._purge_cptrs()
-        rxd.initializer.has_initialized = False
-        rxd.rxd.free_conc_ptrs()
-        rxd.rxd.free_curr_ptrs()
-        rxd.rxd.rxd_include_node_flux1D(0, None, None, None)
-        rxd.species._has_1d = False
-        rxd.species._has_3d = False
-        rxd.rxd._zero_volume_indices = np.ndarray(0, dtype=np.int_)
-        rxd.set_solve_type(dimension=1)
-        # clear reactions in case next sim does not use rxd
-        rxd.rxd.clear_rates()
-
-        for obj in rxd.__dict__:
-            sim.clearObj(obj)
-
-        # except:
-        #    pass
-
     if hasattr(sim, 'net'):
+        if hasattr(sim.net, 'rxd'): # check that 'net' exists before checking sim.net.rxd
+            sim.clearObj(sim.net.rxd)
+            # clean rxd components
+            if 'rxd' not in globals():
+                try:
+                    from neuron import crxd as rxd
+                except:
+                    pass
+            # try:
+            for r in rxd.rxd._all_reactions[:]:
+                if r():
+                    rxd.rxd._unregister_reaction(r)
+
+            for s in rxd.species._all_species:
+                if s():
+                    s().__del__()
+
+            rxd.region._all_regions = []
+            rxd.region._region_count = 0
+            rxd.region._c_region_lookup = None
+            rxd.species._species_counts = 0
+            rxd.section1d._purge_cptrs()
+            rxd.initializer.has_initialized = False
+            rxd.rxd.free_conc_ptrs()
+            rxd.rxd.free_curr_ptrs()
+            rxd.rxd.rxd_include_node_flux1D(0, None, None, None)
+            rxd.species._has_1d = False
+            rxd.species._has_3d = False
+            rxd.rxd._zero_volume_indices = np.ndarray(0, dtype=np.int_)
+            rxd.set_solve_type(dimension=1)
+            # clear reactions in case next sim does not use rxd
+            rxd.rxd.clear_rates()
+
+            for obj in rxd.__dict__:
+                sim.clearObj(obj)
         del sim.net
-
     import gc
-
     gc.collect()
 
+    sim.pc.barrier()
+    sim.pc.gid_clear()  # clear previous gid settings
+
+def close(clear=True):
+    """
+    Function to close simulation
+
+    """
+    from .. import sim
+    if clear:
+        clearAll()
+    else:
+        sim.pc.barrier()
+    sys.exit()
+
 
 def checkConditions(conditions, against, cellGid=None):