diff --git a/mrdna/model/nbPot.py b/mrdna/model/nbPot.py
index e744f82f9b9e0f5eb46db1d90d641dae082300bb..e730a4516e65be7035d3626b20e48edf77967319 100644
--- a/mrdna/model/nbPot.py
+++ b/mrdna/model/nbPot.py
@@ -9,9 +9,6 @@ from mrdna.config import CACHE_DIR
 
 ## TODO: clean this file up a bit, make things protected
 
-_package_cache_dir = get_resource_path("nbPot_cache")
-_package_cache_dir.mkdir(exist_ok=True)
-
 _dxParam = -1
 
 def _maxForce(x,u0,maxForce=40):
diff --git a/mrdna/simulate.py b/mrdna/simulate.py
index df93fc461bf72bb62524d786a98e2b1dba7f8bfb..7494f9b0ed51a28ea77223befe075845006b12b5 100644
--- a/mrdna/simulate.py
+++ b/mrdna/simulate.py
@@ -21,7 +21,9 @@ def multiresolution_simulation( model, output_name,
 
     ## Round steps up to nearest multiple of output_period, plus 1
     coarse_steps = ((coarse_steps+coarse_output_period-1)//coarse_output_period)*coarse_output_period+1
-    fine_steps = ((fine_steps+fine_output_period-1)//fine_output_period)*fine_output_period+1
+    fine_steps = ((fine_steps+fine_output_period-1)//fine_output_period)
+    if fine_steps == 1: fine_steps += 1
+    fine_steps = fine_steps*fine_output_period+1 
 
     ret = None
     d_orig = os.getcwd()
@@ -74,10 +76,16 @@ def multiresolution_simulation( model, output_name,
         model._clear_beads()
         model._generate_bead_model( 1, 1, local_twist=True, escapable_twist=True )
         model.simulate( output_name = output_prefix, num_steps=fine_steps, **simargs )
-        try:
-            coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.dcd' % full_output_prefix, rmsdThreshold=1)
-        except:
-            coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix, rmsdThreshold=1)
+        if fine_steps <= 2*fine_output_period+1:
+            try:
+                coordinates = readArbdCoords('%s.restart' % full_output_prefix)
+            except:
+                coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
+        else:
+            try:
+                coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.dcd' % full_output_prefix, rmsdThreshold=1)
+            except:
+                coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix, rmsdThreshold=1)
 
 
         """ Freeze twist """
@@ -87,10 +95,16 @@ def multiresolution_simulation( model, output_name,
         model._clear_beads()
         model._generate_bead_model( 1, 1, local_twist=True, escapable_twist=False )
         model.simulate( output_name = output_prefix, num_steps=fine_steps, **simargs )
-        try:
-            coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.dcd' % full_output_prefix )
-        except:
-            coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix )
+        if fine_steps <= 2*fine_output_period+1:
+            try:
+                coordinates = readArbdCoords('%s.restart' % full_output_prefix)
+            except:
+                coordinates = readArbdCoords('%s.0.restart' % full_output_prefix)
+        else:
+            try:
+                coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.dcd' % full_output_prefix, rmsdThreshold=1)
+            except:
+                coordinates = readAvgArbdCoords('%s.psf' % output_prefix,'%s.pdb' % output_prefix, '%s.0.dcd' % full_output_prefix, rmsdThreshold=1)
 
 
         """ Atomic simulation """