diff --git a/CHANGES.rst b/CHANGES.rst
index 60dee6c..c304ce7 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -1,6 +1,6 @@
-=========
-Changelog
-=========
+=============
+Release Notes
+=============
 
 0.3.0 (2025-11-03)
 ==================
diff --git a/scheil/simulate.py b/scheil/simulate.py
index e3919aa..f143e1a 100644
--- a/scheil/simulate.py
+++ b/scheil/simulate.py
@@ -51,7 +51,8 @@ def _update_points(wks, points_dict, dof_dict, local_pdens=0, verbose=False):
             if verbose:
                 print(f'Adding points to {phase_name}. ', end='')
             dof = dof_dict[phase_name]
-            eq_pts = np.asarray(compset.dof[:sum(dof)]).reshape(1, -1)
+            num_statevars = compset.phase_record.num_statevars
+            eq_pts = np.asarray(compset.dof[num_statevars:num_statevars + sum(dof)]).reshape(1, -1)
             if local_pdens > 0:
                 points_dict[phase_name] = np.concatenate([pts, local_sample(eq_pts, dof, pdens=local_pdens)], axis=0)
             else:
diff --git a/tests/test_update_points.py b/tests/test_update_points.py
new file mode 100644
index 0000000..93682cb
--- /dev/null
+++ b/tests/test_update_points.py
@@ -0,0 +1,44 @@
+import numpy as np
+from pycalphad import Workspace, variables as v
+from scheil.simulate import _update_points
+from fixtures import select_database, load_database
+
+
+@select_database("alzn_mey.tdb")
+def test_update_points_appends_equilibrium_site_fractions(load_database):
+    """The point appended by _update_points must be the equilibrium site fractions, not the state variables."""
+    dbf = load_database()
+    wks = Workspace(dbf, ['AL', 'ZN', 'VA'], ['LIQUID', 'FCC_A1', 'HCP_A3'], {v.T: 700, v.P: 101325, v.N: 1, v.X('ZN'): 0.3})
+    dof_dict = {ph: list(map(len, mod.constituents)) for ph, mod in wks.models.items()}
+    points_dict = {ph: np.empty((0, sum(dof))) for ph, dof in dof_dict.items()}
+
+    # Single-phase FCC_A1 region for this condition
+    assert [cs.phase_record.phase_name for cs in wks.get_composition_sets()] == ['FCC_A1']
+
+    _update_points(wks, points_dict, dof_dict, local_pdens=0)
+
+    assert points_dict['FCC_A1'].shape == (1, 2)
+    expected = [wks.get(v.Y('FCC_A1', 0, 'AL')), wks.get(v.Y('FCC_A1', 0, 'ZN'))]
+    np.testing.assert_allclose(points_dict['FCC_A1'][0], expected, atol=1e-8)
+    # Other phases were not stable, so nothing should be appended
+    assert points_dict['LIQUID'].shape == (0, 2)
+    assert points_dict['HCP_A3'].shape == (0, 2)
+
+
+@select_database("alzn_mey.tdb")
+def test_update_points_local_samples_are_valid_site_fractions(load_database):
+    """Locally sampled points must be valid site fractions and include the equilibrium point itself."""
+    dbf = load_database()
+    wks = Workspace(dbf, ['AL', 'ZN', 'VA'], ['LIQUID', 'FCC_A1', 'HCP_A3'], {v.T: 700, v.P: 101325, v.N: 1, v.X('ZN'): 0.3})
+    dof_dict = {ph: list(map(len, mod.constituents)) for ph, mod in wks.models.items()}
+    points_dict = {ph: np.empty((0, sum(dof))) for ph, dof in dof_dict.items()}
+
+
+    _update_points(wks, points_dict, dof_dict, local_pdens=20)
+
+    pts = points_dict['FCC_A1']
+    assert pts.shape == (21, 2)  # 20 local samples + the equilibrium point
+    assert np.all(pts >= 0) and np.all(pts <= 1)
+    np.testing.assert_allclose(pts.sum(axis=1), 1.0)
+    expected = [wks.get(v.Y('FCC_A1', 0, 'AL')), wks.get(v.Y('FCC_A1', 0, 'ZN'))]
+    assert np.any(np.all(np.isclose(pts, expected, atol=1e-8), axis=1))