diff --git a/atomistics/workflows/langevin/__init__.py b/atomistics/workflows/langevin/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/atomistics/workflows/langevin/workflow.py b/atomistics/workflows/langevin/workflow.py
new file mode 100644
index 00000000..9302b57c
--- /dev/null
+++ b/atomistics/workflows/langevin/workflow.py
@@ -0,0 +1,137 @@
+import numpy as np
+from scipy.constants import physical_constants
+
+from atomistics.workflows.shared.workflow import Workflow
+
+
+KB = physical_constants["Boltzmann constant in eV/K"][0]
+EV_TO_U_ANGSQ_PER_FSSQ = physical_constants["Faraday constant"][0] / 10**7
+U_ANGSQ_PER_FSSQ_TO_EV = 1.0 / EV_TO_U_ANGSQ_PER_FSSQ
+
+
+def langevin_delta_v(
+    temperature, time_step, masses, velocities, damping_timescale=None
+):
+    """
+    Velocity changes due to the Langevin thermostat.
+    Args:
+        temperature (float): The target temperature in K.
+        time_step (float): The MD time step in fs.
+        masses (numpy.ndarray): Per-atom masses in u with a shape (N_atoms, 1).
+        damping_timescale (float): The characteristic timescale of the thermostat in fs.
+        velocities (numpy.ndarray): Per-atom velocities in angstrom/fs.
+    Returns:
+        (numpy.ndarray): Per atom accelerations to use for changing velocities.
+    """
+    if damping_timescale is not None:
+        drag = -0.5 * time_step * velocities / damping_timescale
+        noise = np.sqrt(
+            EV_TO_U_ANGSQ_PER_FSSQ
+            * KB
+            * temperature
+            * time_step
+            / (masses * damping_timescale)
+        ) * np.random.randn(*velocities.shape)
+        noise -= np.mean(noise, axis=0)
+        return drag + noise
+    else:
+        return 0.0
+
+
+def convert_to_acceleration(forces, masses):
+    return forces * EV_TO_U_ANGSQ_PER_FSSQ / masses
+
+
+def get_initial_velocities(temperature, masses, overheat_fraction=2.0):
+    vel_scale = np.sqrt(EV_TO_U_ANGSQ_PER_FSSQ * KB * temperature / masses) * np.sqrt(
+        overheat_fraction
+    )
+    vel_dir = np.random.randn(len(masses), 3)
+    velocities = vel_scale * vel_dir
+    velocities -= np.mean(velocities, axis=0)
+    return velocities
+
+
+def get_first_half_step(forces, masses, time_step, velocities):
+    acceleration = convert_to_acceleration(forces, masses)
+    return velocities + 0.5 * acceleration * time_step
+
+
+class LangevinWorkflow(Workflow):
+    def __init__(
+        self,
+        structure,
+        temperature=1000.0,
+        overheat_fraction=2.0,
+        damping_timescale=100.0,
+        time_step=1,
+    ):
+        self.structure = structure
+        self.temperature = temperature
+        self.overheat_fraction = overheat_fraction
+        self.damping_timescale = damping_timescale
+        self.time_step = time_step
+        self.masses = np.array([a.mass for a in self.structure[:]])[:, np.newaxis]
+        self.positions = self.structure.positions
+        self.velocities = get_initial_velocities(
+            temperature=self.temperature,
+            masses=self.masses,
+            overheat_fraction=self.overheat_fraction,
+        )
+        self.gamma = self.masses / self.damping_timescale
+        self.forces = None
+
+    def generate_structures(self):
+        """
+
+        Returns:
+            (dict)
+        """
+        if self.forces is not None:
+            # first half step
+            vel_half = get_first_half_step(
+                forces=self.forces,
+                masses=self.masses,
+                time_step=self.time_step,
+                velocities=self.velocities,
+            )
+
+            # damping
+            vel_half += langevin_delta_v(
+                temperature=self.temperature,
+                time_step=self.time_step,
+                masses=self.masses,
+                damping_timescale=self.damping_timescale,
+                velocities=self.velocities,
+            )
+
+            # postion update
+            pos_step = self.positions + vel_half * self.time_step
+            structure = self.structure.copy()
+            structure.positions = pos_step
+        else:
+            structure = self.structure
+        return {"calc_forces": {0: structure}, "calc_energy": {0: structure}}
+
+    def analyse_structures(self, output_dict):
+        self.forces, eng_pot = output_dict["forces"][0], output_dict["energy"][0]
+
+        # second half step
+        acceleration = convert_to_acceleration(forces=self.forces, masses=self.masses)
+        vel_step = self.velocities + 0.5 * acceleration * self.time_step
+
+        # damping
+        vel_step += langevin_delta_v(
+            temperature=self.temperature,
+            time_step=self.time_step,
+            masses=self.masses,
+            damping_timescale=self.damping_timescale,
+            velocities=self.velocities,
+        )
+
+        # kinetic energy
+        kinetic_energy = (
+            0.5 * np.sum(self.masses * vel_step * vel_step) / EV_TO_U_ANGSQ_PER_FSSQ
+        )
+        self.velocities = vel_step.copy()
+        return eng_pot, kinetic_energy
diff --git a/tests/test_langevin_lammps.py b/tests/test_langevin_lammps.py
new file mode 100644
index 00000000..cb908cf1
--- /dev/null
+++ b/tests/test_langevin_lammps.py
@@ -0,0 +1,65 @@
+import os
+
+from ase.build import bulk
+import numpy as np
+import unittest
+
+from atomistics.workflows.langevin.workflow import LangevinWorkflow
+
+
+try:
+    from atomistics.calculators.lammps import (
+        evaluate_with_lammps_library, get_potential_dataframe, LammpsASELibrary
+    )
+
+    skip_lammps_test = False
+except ImportError:
+    skip_lammps_test = True
+
+
+@unittest.skipIf(
+    skip_lammps_test, "LAMMPS is not installed, so the LAMMPS tests are skipped."
+)
+class TestLangevin(unittest.TestCase):
+    def test_langevin(self):
+        steps = 300
+        potential = '1999--Mishin-Y--Al--LAMMPS--ipr1'
+        resource_path = os.path.join(os.path.dirname(__file__), "static", "lammps")
+        structure = bulk("Al", cubic=True).repeat([2, 2, 2])
+        df_pot = get_potential_dataframe(
+            structure=structure,
+            resource_path=resource_path
+        )
+        df_pot_selected = df_pot[df_pot.Name == potential].iloc[0]
+        workflow = LangevinWorkflow(
+            structure=structure,
+            temperature=1000.0,
+            overheat_fraction=2.0,
+            damping_timescale=100.0,
+            time_step=1,
+        )
+        lmp = LammpsASELibrary(
+            working_directory=None,
+            cores=1,
+            comm=None,
+            logger=None,
+            log_file=None,
+            library=None,
+            diable_log_file=True,
+        )
+        eng_pot_lst, eng_kin_lst = [], []
+        for i in range(steps):
+            task_dict = workflow.generate_structures()
+            result_dict = evaluate_with_lammps_library(
+                task_dict=task_dict,
+                potential_dataframe=df_pot_selected,
+                lmp=lmp,
+            )
+            eng_pot, eng_kin = workflow.analyse_structures(output_dict=result_dict)
+            eng_pot_lst.append(eng_pot)
+            eng_kin_lst.append(eng_kin)
+        lmp.close()
+        eng_tot_lst = np.array(eng_pot_lst) + np.array(eng_kin_lst)
+        eng_tot_mean = np.mean(eng_tot_lst[200:])
+        self.assertTrue(-105 < eng_tot_mean)
+        self.assertTrue(eng_tot_mean < -103)