diff --git a/src/atomate2/vasp/flows/core.py b/src/atomate2/vasp/flows/core.py
index d04c58aba6..3dbaec8a7f 100644
--- a/src/atomate2/vasp/flows/core.py
+++ b/src/atomate2/vasp/flows/core.py
@@ -16,7 +16,10 @@
     RelaxMaker,
     StaticMaker,
 )
-from atomate2.vasp.sets.core import HSEBSSetGenerator, NonSCFSetGenerator
+from atomate2.vasp.sets.core import (
+    HSEBSSetGenerator,
+    NonSCFSetGenerator,
+)
 
 if TYPE_CHECKING:
     from pathlib import Path
diff --git a/src/atomate2/vasp/flows/md.py b/src/atomate2/vasp/flows/md.py
new file mode 100644
index 0000000000..ad1a21c29b
--- /dev/null
+++ b/src/atomate2/vasp/flows/md.py
@@ -0,0 +1,166 @@
+"""Flows for running molecular dynamics simulations."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass, field
+from typing import TYPE_CHECKING
+
+from jobflow import Flow, Maker, OutputReference
+
+from atomate2.vasp.jobs.md import MDMaker, md_output
+from atomate2.vasp.sets.core import MDSetGenerator
+
+if TYPE_CHECKING:
+    from pathlib import Path
+
+    from pymatgen.core import Structure
+
+    from atomate2.vasp.jobs.base import BaseVaspMaker
+
+
+@dataclass
+class MultiMDMaker(Maker):
+    """
+    Maker to perform an MD run split in several steps.
+
+    Parameters
+    ----------
+    name : str
+        Name of the flows produced by this maker.
+    md_makers : .BaseVaspMaker
+        Maker to use to generate the first relaxation.
+    """
+
+    name: str = "multi md"
+    md_makers: list[BaseVaspMaker] = field(default_factory=lambda: [MDMaker()])
+
+    def make(
+        self,
+        structure: Structure,
+        prev_dir: str | Path | None = None,
+        prev_traj_ids: list[str] | None = None,
+    ) -> Flow:
+        """
+        Create a flow with several chained MD runs.
+
+        Parameters
+        ----------
+        structure : .Structure
+            A pymatgen structure object.
+        prev_dir : str or Path or None
+            A previous VASP calculation directory to copy output files from.
+        prev_traj_ids: a list of ids of job identifying previous steps of the
+            MD trajectory.
+
+        Returns
+        -------
+        Flow
+            A flow containing n_runs MD calculations.
+        """
+        md_job = None
+        md_jobs = []
+        for i, maker in enumerate(self.md_makers, 1):
+            if md_job is None:
+                md_structure = structure
+                md_prev_dir = prev_dir
+            else:
+                md_structure = md_job.output.structure
+                md_prev_dir = md_job.output.dir_name
+            md_job = maker.make(md_structure, prev_dir=md_prev_dir)
+            md_job.name += f" {i}"
+            md_jobs.append(md_job)
+
+        output_job = md_output(
+            structure=md_jobs[-1].output.structure,
+            vasp_dir=md_jobs[-1].output.dir_name,
+            traj_ids=[j.uuid for j in md_jobs],
+            prev_traj_ids=prev_traj_ids,
+        )
+        output_job.name = "molecular dynamics output"
+
+        md_jobs.append(output_job)
+
+        return Flow(md_jobs, output_job.output, name=self.name)
+
+    def restart_from_uuid(self, md_ref: str | OutputReference) -> Flow:
+        """
+        Create a flow from the output reference of another MultiMDMaker.
+
+        The last output will be used as the starting point and the reference to
+        all the previous steps will be included in the final document.
+
+        Parameters
+        ----------
+        md_ref: str or OutputReference
+            The reference to the output of another MultiMDMaker
+
+        Returns
+        -------
+            A flow containing n_runs MD calculations.
+        """
+        if isinstance(md_ref, str):
+            md_ref = OutputReference(md_ref)
+
+        return self.make(
+            structure=md_ref.structure,
+            prev_dir=md_ref.vasp_dir,
+            prev_traj_ids=md_ref.full_traj_ids,
+        )
+
+    @classmethod
+    def from_parameters(
+        cls,
+        nsteps: int,
+        time_step: float,
+        n_runs: int,
+        ensemble: str,
+        start_temp: float,
+        end_temp: float | None = None,
+        **kwargs,
+    ) -> MultiMDMaker:
+        """
+        Create an instance of the Maker based on the standard parameters.
+
+        Set values in the Flow maker, the Job Maker and the VaspInputGenerator,
+        using them to create the final instance of the Maker.
+
+        Parameters
+        ----------
+        nsteps: int
+            Number of time steps for simulations. The VASP `NSW` parameter.
+        time_step: float
+            The time step (in femtosecond) for the simulation. The VASP
+            `POTIM` parameter.
+        n_runs : int
+            Number of MD runs in the flow.
+        ensemble: str
+            Molecular dynamics ensemble to run. Options include `nvt`, `nve`, and `npt`.
+        start_temp: float
+            Starting temperature. The VASP `TEBEG` parameter.
+        end_temp: float or None
+            Final temperature. The VASP `TEEND` parameter. If None the same
+            as start_temp.
+        kwargs:
+            Other parameters passed
+
+        Returns
+        -------
+            A MultiMDMaker
+        """
+        if end_temp is None:
+            end_temp = start_temp
+        md_makers = []
+        start_temp_i = start_temp
+        increment = (end_temp - start_temp) / n_runs
+        for _ in range(n_runs):
+            end_temp_i = start_temp_i + increment
+            generator = MDSetGenerator(
+                nsteps=nsteps,
+                time_step=time_step,
+                ensemble=ensemble,
+                start_temp=start_temp_i,
+                end_temp=end_temp_i,
+            )
+            md_makers.append(MDMaker(input_set_generator=generator))
+            start_temp_i = end_temp_i
+        return cls(md_makers=md_makers, **kwargs)
diff --git a/src/atomate2/vasp/jobs/core.py b/src/atomate2/vasp/jobs/core.py
index 87b38be08a..a89543588a 100644
--- a/src/atomate2/vasp/jobs/core.py
+++ b/src/atomate2/vasp/jobs/core.py
@@ -6,15 +6,6 @@
 from dataclasses import dataclass, field
 from typing import TYPE_CHECKING
 
-from custodian.vasp.handlers import (
-    FrozenJobErrorHandler,
-    IncorrectSmearingHandler,
-    LargeSigmaHandler,
-    MeshSymmetryErrorHandler,
-    PositiveEnergyErrorHandler,
-    StdErrHandler,
-    VaspErrorHandler,
-)
 from pymatgen.alchemy.materials import TransformedStructure
 from pymatgen.alchemy.transmuters import StandardTransmuter
 
@@ -25,7 +16,6 @@
     HSERelaxSetGenerator,
     HSEStaticSetGenerator,
     HSETightRelaxSetGenerator,
-    MDSetGenerator,
     NonSCFSetGenerator,
     RelaxSetGenerator,
     StaticSetGenerator,
@@ -511,59 +501,3 @@ def make(
         self.write_additional_data.setdefault("transformations:json", tjson)
 
         return super().make.original(self, structure, prev_dir)
-
-
-@dataclass
-class MDMaker(BaseVaspMaker):
-    """
-    Maker to create VASP molecular dynamics jobs.
-
-    Parameters
-    ----------
-    name : str
-        The job name.
-    input_set_generator : .VaspInputSetGenerator
-        A generator used to make the input set.
-    write_input_set_kwargs : dict
-        Keyword arguments that will get passed to :obj:`.write_vasp_input_set`.
-    copy_vasp_kwargs : dict
-        Keyword arguments that will get passed to :obj:`.copy_vasp_outputs`.
-    run_vasp_kwargs : dict
-        Keyword arguments that will get passed to :obj:`.run_vasp`.
-    task_document_kwargs : dict
-        Keyword arguments that will get passed to :obj:`.TaskDoc.from_directory`.
-    stop_children_kwargs : dict
-        Keyword arguments that will get passed to :obj:`.should_stop_children`.
-    write_additional_data : dict
-        Additional data to write to the current directory. Given as a dict of
-        {filename: data}. Note that if using FireWorks, dictionary keys cannot contain
-        the "." character which is typically used to denote file extensions. To avoid
-        this, use the ":" character, which will automatically be converted to ".". E.g.
-        ``{"my_file:txt": "contents of the file"}``.
-    """
-
-    name: str = "molecular dynamics"
-
-    input_set_generator: VaspInputGenerator = field(default_factory=MDSetGenerator)
-
-    # Explicitly pass the handlers to not use the default ones. Some default handlers
-    # such as PotimErrorHandler do not apply to MD runs.
-    run_vasp_kwargs: dict = field(
-        default_factory=lambda: {
-            "handlers": (
-                VaspErrorHandler(),
-                MeshSymmetryErrorHandler(),
-                PositiveEnergyErrorHandler(),
-                FrozenJobErrorHandler(),
-                StdErrHandler(),
-                LargeSigmaHandler(),
-                IncorrectSmearingHandler(),
-            )
-        }
-    )
-
-    # Store ionic steps info in a pymatgen Trajectory object instead of in the output
-    # document.
-    task_document_kwargs: dict = field(
-        default_factory=lambda: {"store_trajectory": True}
-    )
diff --git a/src/atomate2/vasp/jobs/md.py b/src/atomate2/vasp/jobs/md.py
new file mode 100644
index 0000000000..86e4d102e4
--- /dev/null
+++ b/src/atomate2/vasp/jobs/md.py
@@ -0,0 +1,128 @@
+"""Module defining molecular dynamics jobs."""
+
+from __future__ import annotations
+
+import logging
+from dataclasses import dataclass, field
+from typing import TYPE_CHECKING
+
+from custodian.vasp.handlers import (
+    FrozenJobErrorHandler,
+    IncorrectSmearingHandler,
+    LargeSigmaHandler,
+    MeshSymmetryErrorHandler,
+    PositiveEnergyErrorHandler,
+    StdErrHandler,
+    VaspErrorHandler,
+)
+from emmet.core.vasp.calculation import StoreTrajectoryOption
+from jobflow import Response, job
+
+from atomate2.vasp.jobs.base import BaseVaspMaker
+from atomate2.vasp.schemas.md import MultiMDOutput
+from atomate2.vasp.sets.core import MDSetGenerator
+
+if TYPE_CHECKING:
+    from pathlib import Path
+
+    from pymatgen.core import Structure
+
+    from atomate2.vasp.sets.base import VaspInputGenerator
+
+
+logger = logging.getLogger(__name__)
+
+
+@dataclass
+class MDMaker(BaseVaspMaker):
+    """
+    Maker to create VASP molecular dynamics jobs.
+
+    Parameters
+    ----------
+    name : str
+        The job name.
+    input_set_generator : .VaspInputSetGenerator
+        A generator used to make the input set.
+    write_input_set_kwargs : dict
+        Keyword arguments that will get passed to :obj:`.write_vasp_input_set`.
+    copy_vasp_kwargs : dict
+        Keyword arguments that will get passed to :obj:`.copy_vasp_outputs`.
+    run_vasp_kwargs : dict
+        Keyword arguments that will get passed to :obj:`.run_vasp`.
+    task_document_kwargs : dict
+        Keyword arguments that will get passed to :obj:`.TaskDoc.from_directory`.
+    stop_children_kwargs : dict
+        Keyword arguments that will get passed to :obj:`.should_stop_children`.
+    write_additional_data : dict
+        Additional data to write to the current directory. Given as a dict of
+        {filename: data}. Note that if using FireWorks, dictionary keys cannot contain
+        the "." character which is typically used to denote file extensions. To avoid
+        this, use the ":" character, which will automatically be converted to ".". E.g.
+        ``{"my_file:txt": "contents of the file"}``.
+    """
+
+    name: str = "molecular dynamics"
+
+    input_set_generator: VaspInputGenerator = field(default_factory=MDSetGenerator)
+
+    # Explicitly pass the handlers to not use the default ones. Some default handlers
+    # such as PotimErrorHandler do not apply to MD runs.
+    run_vasp_kwargs: dict = field(
+        default_factory=lambda: {
+            "handlers": (
+                VaspErrorHandler(),
+                MeshSymmetryErrorHandler(),
+                PositiveEnergyErrorHandler(),
+                FrozenJobErrorHandler(),
+                StdErrHandler(),
+                LargeSigmaHandler(),
+                IncorrectSmearingHandler(),
+            )
+        }
+    )
+
+    # Store ionic steps info in a pymatgen Trajectory object instead of in the output
+    # document.
+    task_document_kwargs: dict = field(
+        default_factory=lambda: {"store_trajectory": StoreTrajectoryOption.PARTIAL}
+    )
+
+
+@job(output_schema=MultiMDOutput)
+def md_output(
+    structure: Structure,
+    vasp_dir: str | Path,
+    traj_ids: list[str],
+    prev_traj_ids: list[str] | None,
+):
+    """
+    Collect output references of a multistep MD flow.
+
+    Parameters
+    ----------
+    structure: .Structure
+        The final structure to be stored.
+    vasp_dir: str or Path
+        The path to the folder containing the last calculation of a MultiMDMaker.
+    traj_ids: list of str
+        List of the uuids of the jobs that will compose the trajectory.
+    prev_traj_ids: list of str
+        List of the uuids of the jobs coming from previous flow that will be
+        added to the trajectory.
+
+    Returns
+    -------
+    The output dictionary.
+    """
+    full_traj_ids = list(traj_ids)
+    if prev_traj_ids:
+        full_traj_ids = prev_traj_ids + full_traj_ids
+    output = MultiMDOutput.from_structure(
+        structure=structure,
+        meta_structure=structure,
+        vasp_dir=str(vasp_dir),
+        traj_ids=traj_ids,
+        full_traj_ids=full_traj_ids,
+    )
+    return Response(output=output)
diff --git a/src/atomate2/vasp/schemas/md.py b/src/atomate2/vasp/schemas/md.py
new file mode 100644
index 0000000000..c3b6364fb9
--- /dev/null
+++ b/src/atomate2/vasp/schemas/md.py
@@ -0,0 +1,16 @@
+"""Schemas for MD documents."""
+
+from emmet.core.structure import StructureMetadata
+from pydantic import Field
+from pymatgen.core import Structure
+
+
+class MultiMDOutput(StructureMetadata):
+    """Output of a MultiMD Flow."""
+
+    structure: Structure = Field("Final structure of the last step of the flow")
+    vasp_dir: str = Field("Path to the last vasp folder of the flow")
+    traj_ids: list[str] = Field("List of uuids of the MD calculations in the flow")
+    full_traj_ids: list[str] = Field(
+        "List of uuids of the MD calculations in the flow and in previous linked flows"
+    )
diff --git a/src/atomate2/vasp/sets/base.py b/src/atomate2/vasp/sets/base.py
index 397bf2b423..ed0d0ef465 100644
--- a/src/atomate2/vasp/sets/base.py
+++ b/src/atomate2/vasp/sets/base.py
@@ -438,10 +438,20 @@ def get_input_set(
             bandgap=bandgap,
             ispin=ispin,
         )
+        site_properties = structure.site_properties
+        poscar = Poscar(
+            structure,
+            velocities=site_properties.get("velocities"),
+            predictor_corrector=site_properties.get("predictor_corrector"),
+            predictor_corrector_preamble=structure.properties.get(
+                "predictor_corrector_preamble"
+            ),
+            lattice_velocities=structure.properties.get("lattice_velocities"),
+        )
         return VaspInputSet(
             incar=incar,
             kpoints=kpoints,
-            poscar=Poscar(structure),
+            poscar=poscar,
             potcar=self._get_potcar(structure, potcar_spec=potcar_spec),
         )
 
diff --git a/src/atomate2/vasp/sets/core.py b/src/atomate2/vasp/sets/core.py
index a405eb9fd0..bf8996edcf 100644
--- a/src/atomate2/vasp/sets/core.py
+++ b/src/atomate2/vasp/sets/core.py
@@ -823,7 +823,7 @@ class MDSetGenerator(VaspInputGenerator):
     start_temp: float = 300
     end_temp: float = 300
     nsteps: int = 1000
-    time_step: int = 2
+    time_step: float = 2
     auto_ispin: bool = True
 
     def get_incar_updates(
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/inputs/INCAR b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/inputs/INCAR
new file mode 100644
index 0000000000..5385f3acb2
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/inputs/INCAR
@@ -0,0 +1,33 @@
+ALGO = Fast
+EDIFF = 1e-05
+EDIFFG = -0.02
+ENAUG = 1360
+ENCUT = 520
+GGA = Ps
+IBRION = 0
+ISIF = 2
+ISMEAR = 0
+ISPIN = 2
+ISYM = 0
+KBLOCK = 100
+LAECHG = True
+LASPH = True
+LCHARG = False
+LELF = False
+LMIXTAU = True
+LORBIT = 11
+LREAL = False
+LVTOT = True
+LWAVE = False
+MAGMOM = 2*0.6
+MAXMIX = 20
+MDALGO = 2
+NELM = 500
+NELMIN = 4
+NSW = 3
+POTIM = 1
+PREC = Normal
+SIGMA = 0.2
+SMASS = 0
+TEBEG = 300
+TEEND = 300.0
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/inputs/KPOINTS b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/inputs/KPOINTS
new file mode 100644
index 0000000000..093abbe92b
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/inputs/KPOINTS
@@ -0,0 +1,4 @@
+pymatgen with grid density = 100 / number of atoms
+0
+Gamma
+3 3 3
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/inputs/POSCAR b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/inputs/POSCAR
new file mode 100644
index 0000000000..3c72d993a0
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/inputs/POSCAR
@@ -0,0 +1,10 @@
+Si2
+1.0
+   3.3488980000000002    0.0000000000000000    1.9334870000000000
+   1.1162989999999999    3.1573720000000001    1.9334870000000000
+   0.0000000000000000    0.0000000000000000    3.8669750000000001
+Si
+2
+direct
+   0.2500000000000000    0.2500000000000000    0.2500000000000000 Si
+   0.0000000000000000    0.0000000000000000    0.0000000000000000 Si
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/inputs/POTCAR.spec b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/inputs/POTCAR.spec
new file mode 100644
index 0000000000..e267321d2c
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/inputs/POTCAR.spec
@@ -0,0 +1 @@
+Si
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/CONTCAR.gz
new file mode 100644
index 0000000000..6caf0b9685
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/CONTCAR.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/INCAR.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/INCAR.gz
new file mode 100644
index 0000000000..12c4ad3ed3
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/INCAR.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/INCAR.orig.gz
new file mode 100644
index 0000000000..8901363a6b
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/INCAR.orig.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/KPOINTS.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/KPOINTS.gz
new file mode 100644
index 0000000000..d66f7964c9
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/KPOINTS.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/KPOINTS.orig.gz
new file mode 100644
index 0000000000..92b96421da
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/KPOINTS.orig.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/OUTCAR.gz
new file mode 100644
index 0000000000..fa7e41a76c
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/OUTCAR.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/POSCAR.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/POSCAR.gz
new file mode 100644
index 0000000000..b8e86e2924
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/POSCAR.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/POSCAR.orig.gz
new file mode 100644
index 0000000000..2815ed3f58
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/POSCAR.orig.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/POTCAR.spec b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/POTCAR.spec
new file mode 100644
index 0000000000..e267321d2c
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/POTCAR.spec
@@ -0,0 +1 @@
+Si
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/custodian.json.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/custodian.json.gz
new file mode 100644
index 0000000000..1391e055a0
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/custodian.json.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/vasprun.xml.gz
new file mode 100644
index 0000000000..0a59e21721
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_1/outputs/vasprun.xml.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/inputs/INCAR b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/inputs/INCAR
new file mode 100644
index 0000000000..a6eced0c4a
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/inputs/INCAR
@@ -0,0 +1,32 @@
+ALGO = Fast
+EDIFF = 1e-05
+EDIFFG = -0.02
+ENAUG = 1360
+ENCUT = 520
+GGA = Ps
+IBRION = 0
+ISIF = 2
+ISMEAR = -5
+ISPIN = 1
+ISYM = 0
+KBLOCK = 100
+LAECHG = True
+LASPH = True
+LCHARG = False
+LELF = False
+LMIXTAU = True
+LORBIT = 11
+LREAL = False
+LVTOT = True
+LWAVE = False
+MAXMIX = 20
+MDALGO = 2
+NELM = 500
+NELMIN = 4
+NSW = 3
+POTIM = 1
+PREC = Normal
+SIGMA = 0.05
+SMASS = 0
+TEBEG = 300.0
+TEEND = 300.0
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/inputs/KPOINTS b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/inputs/KPOINTS
new file mode 100644
index 0000000000..093abbe92b
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/inputs/KPOINTS
@@ -0,0 +1,4 @@
+pymatgen with grid density = 100 / number of atoms
+0
+Gamma
+3 3 3
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/inputs/POSCAR b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/inputs/POSCAR
new file mode 100644
index 0000000000..d5f7bc1ac2
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/inputs/POSCAR
@@ -0,0 +1,23 @@
+Si2
+1.0
+   3.3488980000000002    0.0000000000000000    1.9334870000000000
+   1.1162989999999999    3.1573720000000001    1.9334870000000000
+   0.0000000000000000    0.0000000000000000    3.8669750000000001
+Si
+2
+direct
+   0.2509246999416777    0.2479629372437218    0.2497071651391614 Si
+  -0.0009246999416777    0.0020370627562782    0.0002928348608386 Si
+
+   0.0004025256800000   -0.0031480439000000   -0.0016045943000000
+  -0.0004025256800000    0.0031480439000000    0.0016045943000000
+
+1
+   1.0000000000000000
+ -0.63509263E-03 -0.84024339E-03 -0.21494156E-03 -0.10750015E-03
+   0.2513772400000000    0.2469658900000000    0.2495644700000000
+   0.9986227600000001    0.0030341084000000    0.0004355327700000
+   0.0000000000000000    0.0000000000000000    0.0000000000000000
+   0.0000000000000000    0.0000000000000000    0.0000000000000000
+   0.0000000000000000    0.0000000000000000    0.0000000000000000
+   0.0000000000000000    0.0000000000000000    0.0000000000000000
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/inputs/POTCAR.spec b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/inputs/POTCAR.spec
new file mode 100644
index 0000000000..e267321d2c
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/inputs/POTCAR.spec
@@ -0,0 +1 @@
+Si
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/CONTCAR.gz
new file mode 100644
index 0000000000..9b64e4cceb
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/CONTCAR.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/INCAR.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/INCAR.gz
new file mode 100644
index 0000000000..cd648b93b4
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/INCAR.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/INCAR.orig.gz
new file mode 100644
index 0000000000..d5276a2199
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/INCAR.orig.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/KPOINTS.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/KPOINTS.gz
new file mode 100644
index 0000000000..1fc7e49e42
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/KPOINTS.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/KPOINTS.orig.gz
new file mode 100644
index 0000000000..65b1aa1275
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/KPOINTS.orig.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/OUTCAR.gz
new file mode 100644
index 0000000000..161b55657b
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/OUTCAR.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/POSCAR.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/POSCAR.gz
new file mode 100644
index 0000000000..08ce6d4fb2
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/POSCAR.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/POSCAR.orig.gz
new file mode 100644
index 0000000000..e8cffe8eb7
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/POSCAR.orig.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/POTCAR.spec b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/POTCAR.spec
new file mode 100644
index 0000000000..e267321d2c
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/POTCAR.spec
@@ -0,0 +1 @@
+Si
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/custodian.json.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/custodian.json.gz
new file mode 100644
index 0000000000..e1b6808bb3
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/custodian.json.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/vasprun.xml.gz
new file mode 100644
index 0000000000..9ed080c665
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_2/outputs/vasprun.xml.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/inputs/INCAR b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/inputs/INCAR
new file mode 100644
index 0000000000..f1fb2a1630
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/inputs/INCAR
@@ -0,0 +1,32 @@
+ALGO = Fast
+EDIFF = 1e-05
+EDIFFG = -0.02
+ENAUG = 1360
+ENCUT = 520
+GGA = Ps
+IBRION = 0
+ISIF = 2
+ISMEAR = -5
+ISPIN = 1
+ISYM = 0
+KBLOCK = 100
+LAECHG = True
+LASPH = True
+LCHARG = False
+LELF = False
+LMIXTAU = True
+LORBIT = 11
+LREAL = False
+LVTOT = True
+LWAVE = False
+MAXMIX = 20
+MDALGO = 2
+NELM = 500
+NELMIN = 4
+NSW = 3
+POTIM = 1
+PREC = Normal
+SIGMA = 0.05
+SMASS = 0
+TEBEG = 300
+TEEND = 300.0
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/inputs/KPOINTS b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/inputs/KPOINTS
new file mode 100644
index 0000000000..093abbe92b
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/inputs/KPOINTS
@@ -0,0 +1,4 @@
+pymatgen with grid density = 100 / number of atoms
+0
+Gamma
+3 3 3
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/inputs/POSCAR b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/inputs/POSCAR
new file mode 100644
index 0000000000..5428183e57
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/inputs/POSCAR
@@ -0,0 +1,23 @@
+Si2
+1.0
+   3.3488980000000002    0.0000000000000000    1.9334870000000000
+   1.1162989999999999    3.1573720000000001    1.9334870000000000
+   0.0000000000000000    0.0000000000000000    3.8669750000000001
+Si
+2
+direct
+   0.2522405717926650    0.2450634335621342    0.2492943829866818 Si
+   0.9977594282073351    0.0049365648378658    0.0007056197833182 Si
+
+   0.0003589612200000   -0.0028134973000000   -0.0014224528000000
+  -0.0003589612200000    0.0028134973000000    0.0014224528000000
+
+1
+   1.0000000000000000
+ -0.10996701E-01 -0.72234638E-02 -0.42099336E-02 -0.52800025E-02
+   0.2526447900000000    0.2441723500000000    0.2491699700000000
+   0.9973552100000000    0.0058276532000000    0.0008300306800000
+   0.0000000000000000    0.0000000000000000    0.0000000000000000
+   0.0000000000000000    0.0000000000000000    0.0000000000000000
+   0.0000000000000000    0.0000000000000000    0.0000000000000000
+   0.0000000000000000    0.0000000000000000    0.0000000000000000
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/inputs/POTCAR.spec b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/inputs/POTCAR.spec
new file mode 100644
index 0000000000..e267321d2c
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/inputs/POTCAR.spec
@@ -0,0 +1 @@
+Si
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/CONTCAR.gz
new file mode 100644
index 0000000000..3fa946cffb
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/CONTCAR.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/INCAR.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/INCAR.gz
new file mode 100644
index 0000000000..da0dc8207a
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/INCAR.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/INCAR.orig.gz
new file mode 100644
index 0000000000..34f8b022d4
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/INCAR.orig.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/KPOINTS.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/KPOINTS.gz
new file mode 100644
index 0000000000..99ba394ea4
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/KPOINTS.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/KPOINTS.orig.gz
new file mode 100644
index 0000000000..bf0816a1f7
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/KPOINTS.orig.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/OUTCAR.gz
new file mode 100644
index 0000000000..6bf737616e
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/OUTCAR.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/POSCAR.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/POSCAR.gz
new file mode 100644
index 0000000000..b57a7ece7b
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/POSCAR.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/POSCAR.orig.gz
new file mode 100644
index 0000000000..4d04e8aca1
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/POSCAR.orig.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/POTCAR.spec b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/POTCAR.spec
new file mode 100644
index 0000000000..e267321d2c
--- /dev/null
+++ b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/POTCAR.spec
@@ -0,0 +1 @@
+Si
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/custodian.json.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/custodian.json.gz
new file mode 100644
index 0000000000..8d4c906334
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/custodian.json.gz differ
diff --git a/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/vasprun.xml.gz
new file mode 100644
index 0000000000..e202d02ca1
Binary files /dev/null and b/tests/test_data/vasp/Si_multi_md/molecular_dynamics_3/outputs/vasprun.xml.gz differ
diff --git a/tests/vasp/flows/test_md.py b/tests/vasp/flows/test_md.py
new file mode 100644
index 0000000000..61478740dc
--- /dev/null
+++ b/tests/vasp/flows/test_md.py
@@ -0,0 +1,70 @@
+from emmet.core.vasp.calculation import VaspObject
+from jobflow import Flow
+
+from atomate2.vasp.flows.md import MultiMDMaker
+from atomate2.vasp.powerups import update_user_kpoints_settings
+from atomate2.vasp.schemas.md import MultiMDOutput
+
+
+def test_multi_md_flow(mock_vasp, clean_dir, si_structure):
+    from emmet.core.tasks import TaskDoc
+    from jobflow import run_locally
+
+    # mapping from job name to directory containing test files
+    ref_paths = {
+        "molecular dynamics 1": "Si_multi_md/molecular_dynamics_1",
+        "molecular dynamics 2": "Si_multi_md/molecular_dynamics_2",
+        "molecular dynamics 3": "Si_multi_md/molecular_dynamics_3",
+    }
+
+    # settings passed to fake_run_vasp; adjust these to check for certain INCAR settings
+    fake_run_vasp_kwargs = {
+        "molecular dynamics 1": {"incar_settings": ["NSW", "ISMEAR"]},
+        "molecular dynamics 2": {"incar_settings": ["NSW", "ISMEAR"]},
+        "molecular dynamics 3": {"incar_settings": ["NSW", "ISMEAR"]},
+    }
+
+    # automatically use fake VASP and write POTCAR.spec during the test
+    mock_vasp(ref_paths, fake_run_vasp_kwargs)
+
+    mdflow1 = MultiMDMaker.from_parameters(
+        nsteps=3, time_step=1, n_runs=2, ensemble="nvt", start_temp=300
+    ).make(si_structure)
+    mdflow2 = MultiMDMaker.from_parameters(
+        nsteps=3, time_step=1, n_runs=1, ensemble="nvt", start_temp=300
+    ).restart_from_uuid(mdflow1.jobs[-1].output)
+    # set the name of the continuation MD Job, otherwise the folders for
+    # files will conflict
+    mdflow2.jobs[0].name = "molecular dynamics 3"
+
+    flow = Flow([mdflow1, mdflow2])
+    flow = update_user_kpoints_settings(
+        flow, {"grid_density": 100}, name_filter="molecular dynamics"
+    )
+
+    # run the flow or job and ensure that it finished running successfully
+    responses = run_locally(flow, create_folders=True, ensure_success=True)
+
+    # validate the outputs
+    output_md_1 = responses[flow.jobs[0].jobs[0].uuid][1].output
+    traj = output_md_1.vasp_objects[VaspObject.TRAJECTORY]
+    assert len(traj.frame_properties) == 3
+    assert isinstance(output_md_1, TaskDoc)
+
+    output_recap_1 = responses[flow.jobs[0].jobs[2].uuid][1].output
+    assert len(output_recap_1.traj_ids) == 2
+    assert len(output_recap_1.full_traj_ids) == 2
+    assert isinstance(output_recap_1, MultiMDOutput)
+
+    output_recap_2 = responses[flow.jobs[1].jobs[1].uuid][1].output
+    assert len(output_recap_2.traj_ids) == 1
+    assert len(output_recap_2.full_traj_ids) == 3
+    assert isinstance(output_recap_1, MultiMDOutput)
+
+
+def test_multi_md_flow_restart_from_uuid():
+    # check that the correct reference is used if a string is passed
+    ref_id = "475bf8ab-06ec-4222-8bad-6f9f3979f2ea"
+    flow = MultiMDMaker().restart_from_uuid(ref_id)
+
+    assert flow.jobs[0].function_kwargs["prev_dir"].uuid == ref_id
diff --git a/tests/vasp/jobs/test_core.py b/tests/vasp/jobs/test_core.py
index 71af761726..633816f069 100644
--- a/tests/vasp/jobs/test_core.py
+++ b/tests/vasp/jobs/test_core.py
@@ -1,8 +1,6 @@
 import jobflow
 import numpy as np
-import pytest
 from emmet.core.tasks import TaskDoc
-from emmet.core.vasp.calculation import IonicStep, VaspObject
 from jobflow import run_locally
 from numpy.testing import assert_allclose
 from pytest import approx
@@ -11,7 +9,6 @@
     DielectricMaker,
     HSERelaxMaker,
     HSEStaticMaker,
-    MDMaker,
     RelaxMaker,
     StaticMaker,
     TransmuterMaker,
@@ -188,48 +185,3 @@ def test_transmuter(mock_vasp, clean_dir, si_structure):
     np.testing.assert_allclose(
         output1.structure.lattice.abc, [3.866974, 3.866975, 7.733949]
     )
-
-
-def test_molecular_dynamics(mock_vasp, clean_dir, si_structure):
-    # mapping from job name to directory containing test files
-    ref_paths = {"molecular dynamics": "Si_molecular_dynamics/molecular_dynamics"}
-
-    # settings passed to fake_run_vasp; adjust these to check for certain INCAR settings
-    fake_run_vasp_kwargs = {
-        "molecular dynamics": {
-            "incar_settings": [
-                "IBRION",
-                "TBEN",
-                "TEND",
-                "NSW",
-                "POTIM",
-                "MDALGO",
-                "ISIF",
-            ]
-        }
-    }
-
-    # automatically use fake VASP and write POTCAR.spec during the test
-    mock_vasp(ref_paths, fake_run_vasp_kwargs)
-
-    # generate job
-    job = MDMaker().make(si_structure)
-    nsw = 3
-    job.maker.input_set_generator.user_incar_settings["NSW"] = nsw
-
-    # run the flow or job and ensure that it finished running successfully
-    responses = run_locally(job, create_folders=True, ensure_success=True)
-
-    # validation on the output
-
-    output1 = responses[job.uuid][1].output
-    assert isinstance(output1, TaskDoc)
-    assert output1.output.energy == pytest.approx(-11.46520398)
-
-    # check ionic steps stored as pymatgen Trajectory
-    assert output1.calcs_reversed[0].output.ionic_steps is None
-    traj = output1.vasp_objects[VaspObject.TRAJECTORY]
-    assert len(traj.frame_properties) == nsw
-    # simply check a frame property can be converted to an IonicStep
-    for frame in traj.frame_properties:
-        IonicStep(**frame)
diff --git a/tests/vasp/jobs/test_md.py b/tests/vasp/jobs/test_md.py
new file mode 100644
index 0000000000..324f360ace
--- /dev/null
+++ b/tests/vasp/jobs/test_md.py
@@ -0,0 +1,53 @@
+import pytest
+from emmet.core.tasks import TaskDoc
+from emmet.core.vasp.calculation import IonicStep, VaspObject
+from jobflow import run_locally
+
+from atomate2.vasp.jobs.md import MDMaker
+
+
+def test_molecular_dynamics(mock_vasp, clean_dir, si_structure):
+    # mapping from job name to directory containing test files
+    ref_paths = {"molecular dynamics": "Si_molecular_dynamics/molecular_dynamics"}
+
+    # settings passed to fake_run_vasp; adjust these to check for certain INCAR settings
+    fake_run_vasp_kwargs = {
+        "molecular dynamics": {
+            "incar_settings": [
+                "IBRION",
+                "TBEN",
+                "TEND",
+                "NSW",
+                "POTIM",
+                "MDALGO",
+                "ISIF",
+            ]
+        }
+    }
+
+    # automatically use fake VASP and write POTCAR.spec during the test
+    mock_vasp(ref_paths, fake_run_vasp_kwargs)
+
+    # generate job
+    job = MDMaker().make(si_structure)
+    nsw = 3
+    job.maker.input_set_generator.user_incar_settings["NSW"] = nsw
+
+    # run the flow or job and ensure that it finished running successfully
+    responses = run_locally(job, create_folders=True, ensure_success=True)
+
+    # validation on the output
+
+    output1 = responses[job.uuid][1].output
+    assert isinstance(output1, TaskDoc)
+    assert output1.output.energy == pytest.approx(-11.46520398)
+
+    # check ionic steps stored as pymatgen Trajectory
+    assert output1.calcs_reversed[0].output.ionic_steps is None
+    traj = output1.vasp_objects[VaspObject.TRAJECTORY]
+    assert len(traj.frame_properties) == nsw
+    # check that a frame property can be converted to an IonicStep
+    energies = [-11.47041923, -11.46905352, -11.46520398]
+    for i, frame in enumerate(traj.frame_properties):
+        ionic_step = IonicStep(**frame)
+        assert ionic_step.e_wo_entrp == pytest.approx(energies[i])