diff --git a/bin/enrgmd b/bin/enrgmd
index 25d9959d3e38ed317be8d781d0d332f3cb8ca396..ed11a1734839de37d9d299be0fa897def973ff74 100755
--- a/bin/enrgmd
+++ b/bin/enrgmd
@@ -78,7 +78,7 @@ http://dx.doi.org/10.1093/nar/gkw155
             os.makedirs('output')
         except FileExistsError:
             ...
-        model.atomic_simulate( output_name = prefix )
+        model.atomic_simulate( output_name = prefix, dry_run=True )
         
     except:
         raise
diff --git a/bin/mrdna b/bin/mrdna
index f5397417fdb9f6670017225ccebd28e23a21996e..06fd281d742c18bfeddea50c35f3d114cc0a1843 100755
--- a/bin/mrdna
+++ b/bin/mrdna
@@ -57,6 +57,10 @@ parser.add_argument('--oxdna-output-period', type=float, default=1e4,
 
 parser.add_argument('--backbone-scale', type=float, default=1.0,
                     help='Factor to scale DNA backbone in atomic model; try 0.25 to avoid clashes for atomistic simulations')
+
+parser.add_argument('--run-enrg-md', action='store_true',
+                    help='Perform the ENRG-MD simulation?')
+
 parser.add_argument('--debug', action='store_true',
                     help='Run through the python debugger?')
 
@@ -128,7 +132,7 @@ def main():
         model._clear_beads()
         model._generate_atomic_model(scale=args.backbone_scale)            
         model.write_atomic_ENM( prefix )
-        model.atomic_simulate( output_name = prefix )
+        model.atomic_simulate( output_name = prefix, dry_run=True )
     elif args.draw_tubes is True:
         output_name = prefix + ".tubes.tcl"
         model.vmd_tube_tcl( output_name )
@@ -140,7 +144,7 @@ def main():
         model._clear_beads()
         model._generate_atomic_model(scale=args.backbone_scale)            
         model.write_atomic_ENM( prefix )
-        model.atomic_simulate( output_name = prefix )
+        model.atomic_simulate( output_name = prefix, dry_run=True )
     else:
         run_args = dict(
             model = model,
@@ -158,7 +162,8 @@ def main():
             fine_steps = int(args.fine_steps),
             backbone_scale = args.backbone_scale,
             oxdna_steps = args.oxdna_steps,
-            oxdna_output_period = args.oxdna_output_period
+            oxdna_output_period = args.oxdna_output_period,
+            run_enrg_md = args.run_enrg_md
         )
 
         simulate( **run_args )
diff --git a/mrdna/model/arbdmodel.py b/mrdna/model/arbdmodel.py
index dd173d2ee5c7d3f9a34fcbd572cfb1a27ab79872..0a358264035998e34dc0ebbd7d07fe1a7cdf11ab 100644
--- a/mrdna/model/arbdmodel.py
+++ b/mrdna/model/arbdmodel.py
@@ -1187,7 +1187,7 @@ if {{$nLast == 0}} {{
 run {num_steps:d}
 """.format(**format_data))
 
-    def atomic_simulate(self, output_name, output_directory='output'):
+    def atomic_simulate(self, output_name, output_directory='output', dry_run = False, namd2=None, log_file=None, num_procs=None, gpu=None):
         if self.cacheUpToDate == False: # TODO: remove cache?
             self._countParticleTypes()
             self._updateParticleOrder()
@@ -1199,6 +1199,48 @@ run {num_steps:d}
         self.write_namd_configuration( output_name, output_directory = output_directory )
         os.sync()
 
+        if not dry_run:
+            if namd2 is None:
+                for path in os.environ["PATH"].split(os.pathsep):
+                    path = path.strip('"')
+                    fname = os.path.join(path, "namd2")
+                    if os.path.isfile(fname) and os.access(fname, os.X_OK):
+                        namd2 = fname
+                        break
+
+            if namd2 is None: raise Exception("NAMD2 was not found")
+
+            if not os.path.exists(namd2):
+                raise Exception("NAMD2 was not found")
+            if not os.path.isfile(namd2):
+                raise Exception("NAMD2 was not found")
+            if not os.access(namd2, os.X_OK):
+                raise Exception("NAMD2 is not executable")
+
+            if not os.path.exists(output_directory):
+                os.makedirs(output_directory)
+            elif not os.path.isdir(output_directory):
+                raise Exception("output_directory '%s' is not a directory!" % output_directory)
+
+            if num_procs is None:
+                import multiprocessing
+                num_procs = max(1,multiprocessing.cpu_count()-1)
+
+            cmd = [namd2, '+p{}'.format(num_procs), "%s.namd" % output_name]
+            cmd = tuple(str(x) for x in cmd)
+
+            print("Running NAMD2 with: %s" % " ".join(cmd))
+            if log_file is None or (hasattr(log_file,'write') and callable(log_file.write)):
+                fd = sys.stdout if log_file is None else log_file
+                process = subprocess.Popen(cmd, stdout=subprocess.PIPE, universal_newlines=True)
+                for line in process.stdout:
+                    fd.write(line)
+                    fd.flush()
+            else:
+                with open(log_file,'w') as fd:
+                    process = subprocess.Popen(cmd, stdout=log_file, universal_newlines=True)
+                    process.communicate()
+
 
     """ oxDNA """
     ## https://dna.physics.ox.ac.uk/index.php/Documentation#Input_file
diff --git a/mrdna/simulate.py b/mrdna/simulate.py
index 348c452e901d0155354857195b25340f8dc40d1a..e8aa7eaebc8469d24f941996820b32eced975f69 100644
--- a/mrdna/simulate.py
+++ b/mrdna/simulate.py
@@ -65,6 +65,7 @@ def multiresolution_simulation( model, output_name,
                                 backbone_scale=1.0,
                                 oxdna_steps = None,
                                 oxdna_output_period = None,
+                                run_enrg_md = False
                             ):
 
     ## Round steps up to nearest multiple of output_period, plus 1
@@ -171,13 +172,14 @@ def multiresolution_simulation( model, output_name,
             # full_output_prefix = "%s/%s" % (output_directory,output_prefix)
 
             model.write_atomic_ENM( output_prefix )
-            model.atomic_simulate( output_name = output_prefix )
-
             try:
                 shutil.copytree( get_resource_path("charmm36.nbfix"), "charmm36.nbfix" )
             except FileExistsError:
                 ...
 
+            model.atomic_simulate( output_name = output_prefix, dry_run = not run_enrg_md )
+
+
 
         ret = directory