diff --git a/matcalc/eos.py b/matcalc/eos.py
new file mode 100644
index 0000000..b2b566f
--- /dev/null
+++ b/matcalc/eos.py
@@ -0,0 +1,71 @@
+"""Calculators for EOS and associated properties."""
+
+from __future__ import annotations
+
+import numpy as np
+from pymatgen.analysis.eos import BirchMurnaghan
+
+from .base import PropCalc
+from .relaxation import RelaxCalc
+
+
+class EOSCalc(PropCalc):
+    """Equation of state calculator."""
+
+    def __init__(
+        self,
+        calculator,
+        relax_structure: bool = True,
+        fmax: float = 0.1,
+        steps: int = 500,
+        n_points: int = 11,
+    ):
+        """
+        Args:
+            calculator: ASE Calculator to use.
+            relax_structure: Whether to first relax the structure. Set to False if structures provided are pre-relaxed
+                with the same calculator.
+            fmax (float): Max force for relaxation (of structure as well as atoms).
+            steps (int): Max number of steps for relaxation.
+            n_points (int): Number of points in which to compute the EOS. Defaults to 11.
+        """
+        self.calculator = calculator
+        self.relax_structure = relax_structure
+        self.n_points = n_points
+        self.fmax = fmax
+        self.steps = steps
+
+    def calc(self, structure):
+        """Fit the Birch-Murnaghan equation of state.
+
+        Args:
+            structure: A Structure object.
+
+        Returns:
+            {
+            "EOS": {
+                "volumes": volumes,
+                "energies": energies
+            },
+            "K": bm.b0,
+        }
+        """
+        if self.relax_structure:
+            relaxer = RelaxCalc(self.calculator, fmax=self.fmax, steps=self.steps)
+            structure = relaxer.calc(structure)["final_structure"]
+
+        volumes, energies = [], []
+        relaxer = RelaxCalc(self.calculator, fmax=self.fmax, steps=self.steps, relax_cell=False)
+        for idx in np.linspace(-0.1, 0.1, self.n_points):
+            structure_strained = structure.copy()
+            structure_strained.apply_strain([idx, idx, idx])
+            result = relaxer.calc(structure_strained)
+            volumes.append(result["final_structure"].volume)
+            energies.append(result["energy"])
+        bm = BirchMurnaghan(volumes=volumes, energies=energies)
+        bm.fit()
+
+        return {
+            "EOS": {"volumes": volumes, "energies": energies},
+            "K (GPa)": bm.b0_GPa,
+        }
diff --git a/matcalc/relaxation.py b/matcalc/relaxation.py
index 235422e..a562313 100644
--- a/matcalc/relaxation.py
+++ b/matcalc/relaxation.py
@@ -92,16 +92,18 @@ def __init__(
         steps: int = 500,
         traj_file: str | None = None,
         interval=1,
+        relax_cell=True,
     ):
         """
         Args:
             calculator: ASE Calculator to use.
             optimizer (str or ase Optimizer): the optimization algorithm.
                 Defaults to "FIRE"
-            fmax (float): total force tolerance for relaxation convergence. fmax is a sum of force and stress forces.
-            steps (int): max number of steps for relaxation.
-            traj_file (str): the trajectory file for saving
-            interval (int): the step interval for saving the trajectories.
+            fmax (float): Total force tolerance for relaxation convergence. fmax is a sum of force and stress forces.
+            steps (int): Max number of steps for relaxation.
+            traj_file (str): The trajectory file for saving
+            interval (int): The step interval for saving the trajectories.
+            relax_cell (bool): Whether to relax the cell.
         """
         self.calculator = calculator
         self.optimizer: Optimizer = OPTIMIZERS[optimizer] if isinstance(optimizer, str) else optimizer
@@ -109,6 +111,7 @@ def __init__(
         self.interval = interval
         self.steps = steps
         self.traj_file = traj_file
+        self.relax_cell = relax_cell
 
     def calc(self, structure) -> dict:
         """
@@ -134,14 +137,16 @@ def calc(self, structure) -> dict:
         stream = io.StringIO()
         with contextlib.redirect_stdout(stream):
             obs = TrajectoryObserver(atoms)
-            atoms = ExpCellFilter(atoms)
+            if self.relax_cell:
+                atoms = ExpCellFilter(atoms)
             optimizer = self.optimizer(atoms)
             optimizer.attach(obs, interval=self.interval)
             optimizer.run(fmax=self.fmax, steps=self.steps)
             if self.traj_file is not None:
                 obs()
                 obs.save(self.traj_file)
-        atoms = atoms.atoms
+        if self.relax_cell:
+            atoms = atoms.atoms
 
         final_structure = ase_adaptor.get_structure(atoms)
         lattice = final_structure.lattice
@@ -155,4 +160,5 @@ def calc(self, structure) -> dict:
             "beta": lattice.beta,
             "gamma": lattice.gamma,
             "volume": lattice.volume,
+            "energy": obs.energies[-1],
         }
diff --git a/tests/test_eos.py b/tests/test_eos.py
new file mode 100644
index 0000000..f5d0239
--- /dev/null
+++ b/tests/test_eos.py
@@ -0,0 +1,20 @@
+"""Tests for PhononCalc class"""
+from __future__ import annotations
+
+import pytest
+
+from matcalc.eos import EOSCalc
+
+
+def test_PhononCalc(Li2O, LiFePO4, M3GNetUPCalc):
+    """Tests for PhononCalc class"""
+    calculator = M3GNetUPCalc
+    # Note that the fmax is probably too high. This is for testing purposes only.
+    pcalc = EOSCalc(calculator)
+    results = pcalc.calc(Li2O)
+
+    assert results["K (GPa)"] == pytest.approx(69.86879801931632, 1e-4)
+
+    results = list(pcalc.calc_many([Li2O, LiFePO4]))
+    assert len(results) == 2
+    assert results[1]["K (GPa)"] == pytest.approx(60.083101510774476, 1e-4)