diff --git a/src/pyiron_lammps/compatibility/file.py b/src/pyiron_lammps/compatibility/file.py
index 709b873..d072d0f 100644
--- a/src/pyiron_lammps/compatibility/file.py
+++ b/src/pyiron_lammps/compatibility/file.py
@@ -1,4 +1,5 @@
 import os
+import shutil
 import subprocess
 from typing import Optional
 
@@ -26,6 +27,9 @@ def lammps_file_interface_function(
     lmp_command: str = "mpiexec -n 1 --oversubscribe lmp_mpi -in lmp.in",
     resource_path: Optional[str] = None,
     input_control_file: Optional[dict] = None,
+    write_restart_file: bool = False,
+    read_restart_file: bool = False,
+    restart_file: str = "restart.out",
 ):
     """
     A single function to execute a LAMMPS calculation based on the LAMMPS job implemented in pyiron
@@ -94,7 +98,12 @@ def lammps_file_interface_function(
 
     lmp_str_lst = []
     atom_type = "atomic"
-    for l in lammps_file_initialization(structure=structure):
+    for l in lammps_file_initialization(
+        structure=structure,
+        units=units,
+        read_restart_file=read_restart_file,
+        restart_file=restart_file,
+    ):
         if l.startswith("units") and "units" in potential_replace:
             lmp_str_lst.append(potential_replace["units"])
         elif l.startswith("atom_style") and "atom_style" in potential_replace:
@@ -136,12 +145,16 @@ def lammps_file_interface_function(
             n_ionic_steps = int(calc_kwargs.pop("n_ionic_steps"))
         else:
             n_ionic_steps = 1
+        if read_restart_file:
+            calc_kwargs["initial_temperature"] = 0.0
         calc_kwargs["units"] = units
         lmp_str_lst += calc_md(**calc_kwargs)
         lmp_str_lst = _modify_input_dict(
             input_control_file=input_control_file,
             lmp_str_lst=lmp_str_lst,
         )
+        if read_restart_file:
+            lmp_str_lst += ["reset_timestep 0"]
         lmp_str_lst += ["run {} ".format(n_ionic_steps)]
     elif calc_mode == "minimize":
         calc_kwargs["units"] = units
@@ -161,6 +174,13 @@ def lammps_file_interface_function(
         raise ValueError(
             f"calc_mode must be one of: static, md or minimize, not {calc_mode}"
         )
+    if read_restart_file:
+        shutil.copyfile(
+            os.path.abspath(restart_file),
+            os.path.join(working_directory, os.path.basename(restart_file)),
+        )
+    if write_restart_file:
+        lmp_str_lst.append(f"write_restart {os.path.basename(restart_file)}")
 
     with open(os.path.join(working_directory, "lmp.in"), "w") as f:
         f.writelines([l + "\n" for l in lmp_str_lst])
@@ -195,15 +215,20 @@ def lammps_file_interface_function(
     return shell, output, False
 
 
-def lammps_file_initialization(structure, dimension=3, units="metal"):
+def lammps_file_initialization(
+    structure, dimension=3, units="metal", read_restart_file=False, restart_file=None
+):
+    init_commands = ["units " + units]
     boundary = " ".join(["p" if coord else "f" for coord in structure.pbc])
-    init_commands = [
-        "units " + units,
-        "dimension " + str(dimension),
-        "boundary " + boundary + "",
-        "atom_style atomic",
-        "read_data lammps.data",
-    ]
+    if read_restart_file:
+        init_commands.append(f"read_restart {os.path.basename(restart_file)}")
+    else:
+        init_commands += [
+            "dimension " + str(dimension),
+            "boundary " + boundary + "",
+            "atom_style atomic",
+            "read_data lammps.data",
+        ]
     return init_commands
 
 
diff --git a/tests/static/restart/restart.out b/tests/static/restart/restart.out
new file mode 100644
index 0000000..e69de29
diff --git a/tests/test_compatibility_file.py b/tests/test_compatibility_file.py
index 93cef2b..c1e4b6b 100644
--- a/tests/test_compatibility_file.py
+++ b/tests/test_compatibility_file.py
@@ -230,6 +230,48 @@ def test_calc_md_nvt(self):
         for line in content_expected:
             self.assertIn(line, content)
 
+    def test_calc_md_nvt_restart(self):
+        calc_kwargs = {"temperature": 500.0, "n_print": 100}
+        shell_output, parsed_output, job_crashed = lammps_file_interface_function(
+            working_directory=self.working_dir,
+            structure=self.structure,
+            potential=self.potential,
+            calc_mode="md",
+            calc_kwargs=calc_kwargs,
+            units=self.units,
+            lmp_command="cp "
+            + str(os.path.join(self.static_path, "compatibility_output"))
+            + "/* .",
+            resource_path=os.path.join(self.static_path, "potential"),
+            restart_file=os.path.join(self.static_path, "restart", "restart.out"),
+            write_restart_file=True,
+            read_restart_file=True,
+        )
+        self.assertFalse(job_crashed)
+        for key in self.keys:
+            self.assertIn(key, parsed_output["generic"])
+        with open(self.working_dir + "/lmp.in", "r") as f:
+            content = f.readlines()
+        content_expected = [
+            "units metal\n",
+            "read_restart restart.out\n",
+            "pair_style eam/alloy\n",
+            "variable dumptime equal 100 \n",
+            "dump 1 all custom ${dumptime} dump.out id type xsu ysu zsu fx fy fz vx vy vz\n",
+            'dump_modify 1 sort id format line "%d %d %20.15g %20.15g %20.15g %20.15g %20.15g %20.15g %20.15g %20.15g %20.15g"\n',
+            "fix ensemble all nvt temp 500.0 500.0 0.1\n",
+            "variable thermotime equal 100 \n",
+            "timestep 0.001\n",
+            "thermo_style custom step temp pe etotal pxx pxy pxz pyy pyz pzz vol\n",
+            "thermo_modify format float %20.15g\n",
+            "thermo ${thermotime}\n",
+            "reset_timestep 0\n",
+            "run 1 \n",
+            "write_restart restart.out\n",
+        ]
+        for line in content_expected:
+            self.assertIn(line, content)
+
     def test_calc_md_nvt_langevin(self):
         calc_kwargs = {"temperature": 500.0, "n_print": 100, "langevin": True}
         shell_output, parsed_output, job_crashed = lammps_file_interface_function(