diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 0000000..ef9d0d5
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,78 @@
+name: CI
+
+on:
+  push:
+    branches: [main]
+  pull_request:
+  workflow_dispatch:
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  build:
+    name: ${{ matrix.os }} / Python ${{ matrix.python-version }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-latest, macos-latest]
+        python-version: ["3.10", "3.12"]
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+
+      - name: Install build dependencies (Linux)
+        if: runner.os == 'Linux'
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y gfortran liblapack-dev libblas-dev pkg-config
+
+      - name: Install build dependencies (macOS)
+        if: runner.os == 'macOS'
+        run: |
+          brew install gcc pkg-config
+          # gfortran ships inside the gcc formula
+          echo "FC=$(brew --prefix gcc)/bin/gfortran" >> "$GITHUB_ENV"
+
+      - name: Build and install
+        run: python -m pip install --upgrade pip && python -m pip install '.[test]'
+
+      - name: Check the console scripts are installed
+        run: |
+          csld_main --help
+          cs-fit --help
+          polaron_main --help
+
+      - name: Run tests
+        run: python -m pytest tests -v
+
+  wheel:
+    name: Build wheel
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+
+      - name: Install build dependencies
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y gfortran liblapack-dev libblas-dev pkg-config
+
+      - name: Build wheel
+        run: |
+          python -m pip install --upgrade pip build
+          python -m build --wheel
+
+      - uses: actions/upload-artifact@v4
+        with:
+          name: csld-wheel
+          path: dist/*.whl
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..1ccd1fd
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,24 @@
+# build artefacts
+build/
+builddir/
+dist/
+*.egg-info/
+.mesonpy-*/
+
+# compiled objects
+*.so
+*.o
+*.a
+*.mod
+*.pyc
+__pycache__/
+
+# cython/f2py generated sources
+compile/c_util/_c_util.cpp
+*module.c
+*-f2pywrappers2.f90
+.f2py_f2cmap
+
+# tooling
+.pytest_cache/
+.venv/
diff --git a/MANIFEST.in b/MANIFEST.in
deleted file mode 100644
index 9948332..0000000
--- a/MANIFEST.in
+++ /dev/null
@@ -1 +0,0 @@
-include compile/c_util/*.h
diff --git a/Makefile b/Makefile
deleted file mode 100644
index 2492cc7..0000000
--- a/Makefile
+++ /dev/null
@@ -1,97 +0,0 @@
-UNAME_S := $(shell uname -s)
-UNAME_N := $(shell uname -n | sed 's/[0-9]*//g')
-ARCH := $(shell uname -p)
-
-
-ifeq ($(UNAME_S),Linux)
-ifeq ($(ARCH),ppc64le)
-LIB =
-F2PY = f2py3
-OPTFLAG = -O
-FFLAGS= -fopenmp -cpp
-CMD = ${F2PY} --fcompiler=gfortran -c  --link-lapack_opt --opt="-cpp -heap-arrays ${OPTFLAG}" ${LIB} --f90flags="${FFLAGS}"
-F90=gfortran
-F77=${F90}
-F77FLAGS = ${FFLAGS}
-else
-# Use the active Python's f2py and a GNU Fortran toolchain by default.
-# These variables can still be overridden, e.g. `make F90=ifort F2PY=f2py3`.
-LIB ?=
-F2PY ?= python3 -m numpy.f2py
-OPTFLAG ?= -O3
-FFLAGS ?= -fopenmp -cpp
-CMD = ${F2PY} -c --link-lapack_opt --opt="${OPTFLAG}" ${LIB} --f90flags="${FFLAGS}"
-F90 ?= gfortran
-F77 = ${F90}
-F77FLAGS ?= ${FFLAGS}
-endif
-endif
-
-ifeq ($(UNAME_S),Darwin)
-LIB = -L${HOME}/lib
-F2PY = f2py
-OPTFLAG = -O
-FFLAGS = -ffree-line-length-none -std=f2008 -fopenmp -cpp
-CMD = ${F2PY} --fcompiler=gfortran -c  --link-lapack_opt --opt="-cpp ${OPTFLAG}" ${LIB} --f90flags="${FFLAGS}"
-F90 = gfortran
-# WARNING: DO NOT use -O* options otherwise code behavior is crazy, seemingly due to the associate statement of fortran 2003
-F77=${F90}
-F77FLAGS =-fopenmp -O -cpp
-endif
-#FFLAGS+= -g -fbacktrace -fbounds-check -ffree-line-length-none -ftrapv
-
-.SUFFIXES:
-
-#BINARIES = bregman.so f_phonon.so f_util.so bcs_driver.so
-BINARIES = bregman.so f_phonon.so f_util.so
-
-BCSSOLVER= cssolve/num_types.f90  cssolve/matrix_sets.f90  cssolve/laplace.f90  cssolve/bcs.f90
-BCSSOLVEROBJ=$(BCSSOLVER:.f90=.o)
-
-
-TET = csld/phonon/dostet.f90 csld/phonon/pdstet.f90 csld/phonon/setk06.f90
-TET_OBJ=$(TET:.f90=.o)
-DMdipole = csld/phonon/DM_dipole_dipole.f csld/phonon/gbox.f csld/phonon/tripl.f csld/phonon/derfc.f csld/phonon/find_Ewald_eta_screened.f csld/phonon/DMstandard.f csld/phonon/lattc.f csld/phonon/cross.f csld/phonon/Fewald.f
-DMdipole_OBJ=$(DMdipole:.f=.o)
-
-.PHONY: all FORCE
-all : $(BINARIES)
-
-$(TET_OBJ) $(DMdipole_OBJ): FORCE
-
-bregman.so: cssolve/bregman.f90
-
-f_phonon.so: csld/phonon/f_phonon.f90 ${TET_OBJ} ${DMdipole_OBJ}
-
-f_util.so: compile/f_util/f_util.f90
-
-clean:
-	rm -f cssolve/*.o cssolve/*.mod cssolve/*.a *.o *.mod $(BINARIES) csld/phonon/*.o
-
-
-%.o: %.f90
-	$(F90) -fPIC -c $(FFLAGS) $< -o $@
-
-%.o: %.f
-	$(F77) -fPIC -c $(F77FLAGS) $< -o $@
-
-cssolve/bcs.a: $(BCSSOLVEROBJ)
-	ar ru $@ $?
-	ranlib $@
-
-cssolve/bcs.so: $(BCSSOLVEROBJ)
-	$(F90) -shared -fPIC $(FFLAGS) -o $@ $^
-
-
-$(BINARIES):
-	rm -f $@
-	${CMD} -m $(basename $@) $^
-	if [ ! -f $@ ]; then mv -f $(basename $@).*.so $@; fi
-
-
-bcs_driver.so: cssolve/bcs_driver.f90 $(BCSSOLVEROBJ)
-	rm -f $@
-	echo "dict(real=dict(dp='double'))" > .f2py_f2cmap
-	${CMD} -m $(basename $@) $^
-	if [ ! -f $@ ]; then mv -f $(basename $@).*.so $@; fi
-	rm .f2py_f2cmap
diff --git a/PKG-INFO b/PKG-INFO
deleted file mode 100644
index 970ede3..0000000
--- a/PKG-INFO
+++ /dev/null
@@ -1,10 +0,0 @@
-Metadata-Version: 1.0
-Name: csld
-Version: 1.0
-Summary: CSLD
-Home-page: https://github.com/LLNL/csld/
-Author: Fei Zhou
-Author-email: fei.fzhou@gmail.com
-License: MIT
-Description: Compressive sensing lattice dynamics
-Platform: any
diff --git a/README.md b/README.md
index 463c64f..b757dc7 100644
--- a/README.md
+++ b/README.md
@@ -4,12 +4,41 @@ A comprehensive package to study lattice dynamics from first-principles. The int
 Getting Started
 -------------------
 
-To install CSLD, make sure you have Python, C++ and Fortran 90 compilers.
-Then:
+CSLD builds C++ and Fortran extensions, so you need a C++ compiler, a Fortran 90
+compiler and BLAS/LAPACK. Everything else is handled by `pip`.
 
     $ git clone https://github.com/LLNL/csld.git
     $ cd csld
-    $ ./setup.py install --user
+    $ pip install .
+
+### With conda
+
+If you do not have compilers or BLAS/LAPACK available, the bundled environment
+provides them:
+
+    $ conda env create -f environment.yml
+    $ conda activate csld
+    $ pip install .
+
+### With a system toolchain
+
+On Debian/Ubuntu, `sudo apt install gfortran liblapack-dev libblas-dev`; on
+macOS, `brew install gcc lapack`. Then `pip install .` as above. Use
+`pip install --user .` to install without admin rights.
+
+### For development
+
+    $ pip install -e . --no-build-isolation
+
+An editable install serves the Python modules straight from the checkout and
+rebuilds the compiled extensions on import, replacing the old
+`setup.py develop`.
+
+After installing, the `csld_main`, `cs-fit`, `polaron_main` and `phonopy-qha`
+commands are on your `PATH`. To check the install:
+
+    $ pip install '.[test]'
+    $ pytest
 
 Documentation
 ----------------
diff --git a/compile/c_util/_c_util.pyx b/compile/c_util/_c_util.pyx
index 560e312..8b06418 100644
--- a/compile/c_util/_c_util.pyx
+++ b/compile/c_util/_c_util.pyx
@@ -1,3 +1,4 @@
+# cython: language_level=3
 """
 C code for csld, including
 fct_trans
diff --git a/compile/c_util/meson.build b/compile/c_util/meson.build
new file mode 100644
index 0000000..84e13b0
--- /dev/null
+++ b/compile/c_util/meson.build
@@ -0,0 +1,20 @@
+# Cython/C++ extension, installed as csld._c_util
+py.extension_module(
+  '_c_util',
+  [
+    '_c_util.pyx',
+    'fct_trans.cpp',
+    'iid-mc.cpp',
+    'intDigits.cpp',
+    'LDCorr.cpp',
+    'LDFFCorr.cpp',
+    'nullspace.cpp',
+    'structure_ordering.cpp',
+  ],
+  override_options: ['cython_language=cpp'],
+  cpp_args: ['-DNDEBUG'],
+  include_directories: [inc_np, include_directories('.')],
+  dependencies: [py_dep],
+  install: true,
+  subdir: 'csld',
+)
diff --git a/compile/f_util/meson.build b/compile/f_util/meson.build
new file mode 100644
index 0000000..968a0b7
--- /dev/null
+++ b/compile/f_util/meson.build
@@ -0,0 +1,17 @@
+# f2py extension, installed as csld.f_util
+f_util_gen = custom_target(
+  'f_utilmodule.c',
+  input: ['f_util.f90'],
+  output: ['f_utilmodule.c', 'f_util-f2pywrappers2.f90'],
+  command: [f2py, '@INPUT@', '-m', 'f_util', '--lower', '--build-dir', '@OUTDIR@'],
+)
+
+py.extension_module(
+  'f_util',
+  ['f_util.f90', f_util_gen, fortranobject_c],
+  include_directories: [inc_np, inc_f2py],
+  dependencies: [py_dep],
+  fortran_args: fortran_args,
+  install: true,
+  subdir: 'csld',
+)
diff --git a/csld/analyzer.py b/csld/analyzer.py
index b876ec9..09fecfc 100644
--- a/csld/analyzer.py
+++ b/csld/analyzer.py
@@ -77,7 +77,7 @@ def get_space_group_symbol(self):
         Returns:
             (str): Spacegroup symbol for structure.
         """
-        return self._space_group_data["international"]
+        return self._space_group_data.international
 
     def get_space_group_number(self):
         """
@@ -86,7 +86,7 @@ def get_space_group_number(self):
         Returns:
             (int): International spacegroup number for structure.
         """
-        return int(self._space_group_data["number"])
+        return int(self._space_group_data.number)
 
     def get_hall(self):
         """
@@ -95,7 +95,7 @@ def get_hall(self):
         Returns:
             (str): Hall symbol
         """
-        return self._space_group_data["hall"]
+        return self._space_group_data.hall
 
     def get_crystal_system(self):
         """
@@ -105,7 +105,7 @@ def get_crystal_system(self):
         Returns:
             (str): Crystal system for structure.
         """
-        n = self._space_group_data["number"]
+        n = self._space_group_data.number
 
         f = lambda i, j: i <= n <= j
         cs = {"triclinic": (1, 2), "monoclinic": (3, 15),
@@ -131,7 +131,7 @@ def get_lattice_type(self):
         Returns:
             (str): Lattice type for structure.
         """
-        n = self._space_group_data["number"]
+        n = self._space_group_data.number
         system = self.get_crystal_system()
         if n in [146, 148, 155, 160, 161, 166, 167]:
             return "rhombohedral"
@@ -142,10 +142,10 @@ def get_lattice_type(self):
 
     def get_symmetry_dataset(self):
         """
-        Returns the symmetry dataset as a dict.
+        Returns the symmetry dataset as returned by spglib.
 
         Returns:
-            (dict): With the following properties:
+            (spglib.SpglibDataset): With the following attributes:
             number: International space group number
             international: International symbol
             hall: Hall symbol
@@ -219,7 +219,7 @@ def get_point_group_symbol(rotations):
         return s[0].strip(), s[1]
 
     def point_group_symbol(self):
-        return SpacegroupAnalyzer.get_point_group_symbol(self._space_group_data["rotations"])
+        return SpacegroupAnalyzer.get_point_group_symbol(self._space_group_data.rotations)
 
     def get_point_group_operations(self):
         """
@@ -254,7 +254,7 @@ def get_symmetrized_structure(self, check_prim=2):
             :class:`SymmetrizedStructure` object.
         """
         return SymmetrizedStructure(self._structure, self.get_spacegroup(),
-                                    self.get_symmetry_dataset()["equivalent_atoms"], self, warn_prim=check_prim)
+                                    self.get_symmetry_dataset().equivalent_atoms, self, warn_prim=check_prim)
 
     def get_refined_structure(self):
         """
diff --git a/csld/basic_lattice_model.py b/csld/basic_lattice_model.py
index 9f3fefe..f45def7 100644
--- a/csld/basic_lattice_model.py
+++ b/csld/basic_lattice_model.py
@@ -14,7 +14,7 @@
 from .atomic_model import AtomicModel
 from .util.mathtool import relativePosition
 from .util.tool import pad_right, matrix2text
-from _c_util import get_structure_ordering
+from csld._c_util import get_structure_ordering
 logger = logging.getLogger(__name__)
 
 
@@ -274,7 +274,7 @@ def translate_cluster_to_supercell(sc, clus):
         coords= clus.frac_coords
         use_compiled_code = True
         if use_compiled_code:
-            from f_util import f_util
+            from csld.f_util import f_util
             allclus= f_util.tr_cl_sc(sc._ijkl.T, sc.sc_mat.T, sc.inv_sc_mat.T, sc.sc_ref.T, sc.prim.frac_coords.T, clus._ijkls_np.T).T-1
         else:
             allclus = [get_structure_ordering((coords+ijk[None,:]).dot(sc.inv_sc_mat), sc.frac_coords,0) for ijk in sc.ijk_ref]
diff --git a/csld/cli/__init__.py b/csld/cli/__init__.py
new file mode 100644
index 0000000..c94a6bc
--- /dev/null
+++ b/csld/cli/__init__.py
@@ -0,0 +1,14 @@
+"""Command line entry points for CSLD.
+
+Each module in this package exposes a ``main()`` function that is registered as a
+console script in ``pyproject.toml``:
+
+===================  ==========================
+console script       module
+===================  ==========================
+``csld_main``        :mod:`csld.cli.csld_main`
+``cs-fit``           :mod:`csld.cli.cs_fit`
+``polaron_main``     :mod:`csld.cli.polaron_main`
+``phonopy-qha``      :mod:`csld.cli.phonopy_qha`
+===================  ==========================
+"""
diff --git a/scripts/cs-fit b/csld/cli/cs_fit.py
similarity index 93%
rename from scripts/cs-fit
rename to csld/cli/cs_fit.py
index d83ba7a..46c15e8 100755
--- a/scripts/cs-fit
+++ b/csld/cli/cs_fit.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+"""``cs-fit`` -- solve a linear model A x = b with compressive sensing."""
 
 import numpy as np
 import sys
@@ -45,7 +45,7 @@ def init_cmdline_settingfile():
 
 
 
-if __name__ == '__main__':
+def main():
     options = init_cmdline_settingfile()
 
     Amat = scipy.io.mmread(options.Amat)
@@ -60,7 +60,8 @@ def init_cmdline_settingfile():
     if options.file == '':
         # simple, oneshot fitting
         print('debug subm=', options.submodel)
-        ibest, solutions = csfit(Amat, fval[:,0], 1, mulist,
+        # csfit() returns (ibest, solutions, rel_err); cf. csld.common_main.fit_data
+        ibest, solutions, rel_err = csfit(Amat, fval[:,0], 1, mulist,
                 method=options.method,
                 maxIter=500,
                 tol=0.0001,
@@ -87,3 +88,7 @@ def init_cmdline_settingfile():
 
     pdfout.close()
 
+
+if __name__ == '__main__':
+    main()
+
diff --git a/scripts/csld_main b/csld/cli/csld_main.py
similarity index 97%
rename from scripts/csld_main
rename to csld/cli/csld_main.py
index 45413d8..ac22ba3 100755
--- a/scripts/csld_main
+++ b/csld/cli/csld_main.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+"""``csld_main`` -- fit interatomic force constants and compute phonons."""
 
 import sys
 import os
@@ -83,7 +83,7 @@ def init_cmdline_settingfile():
     return options, settings
 
 
-if __name__ == '__main__':
+def main():
     logger = logging.getLogger(__name__)
     options, settings = init_cmdline_settingfile()
     logger.setLevel(level=options.log_level)
@@ -102,3 +102,7 @@ def init_cmdline_settingfile():
     if settings.has_section('prediction'):
         predict(model, solutions, settings['prediction'], options.pred_step)
 
+
+if __name__ == '__main__':
+    main()
+
diff --git a/csld/cli/meson.build b/csld/cli/meson.build
new file mode 100644
index 0000000..8aafecf
--- /dev/null
+++ b/csld/cli/meson.build
@@ -0,0 +1,10 @@
+py.install_sources(
+  [
+    '__init__.py',
+    'cs_fit.py',
+    'csld_main.py',
+    'phonopy_qha.py',
+    'polaron_main.py',
+  ],
+  subdir: 'csld/cli',
+)
diff --git a/scripts/phonopy-qha b/csld/cli/phonopy_qha.py
similarity index 93%
rename from scripts/phonopy-qha
rename to csld/cli/phonopy_qha.py
index 0ef5679..2ab2651 100755
--- a/scripts/phonopy-qha
+++ b/csld/cli/phonopy_qha.py
@@ -36,9 +36,17 @@
 
 import numpy as np
 import sys
-from phonopy.units import EVAngstromToGPa
-from phonopy import PhonopyQHA
-from phonopy.file_IO import read_thermal_properties_yaml, read_v_e, read_efe
+
+# phonopy is an optional dependency of csld (`pip install 'csld[qha]'`); defer
+# reporting a missing install until the console script is actually run.
+try:
+    from phonopy.units import EVAngstromToGPa
+    from phonopy import PhonopyQHA
+    from phonopy.file_IO import read_thermal_properties_yaml, read_v_e, read_efe
+except ImportError as exc:  # pragma: no cover - depends on optional extra
+    _phonopy_import_error = exc
+else:
+    _phonopy_import_error = None
 
 
 # simplified file reader assuming plain text format for all data
@@ -109,7 +117,7 @@ def get_options():
     return args
 
 
-def main(args):
+def run(args):
     if args.is_graph_save:
         import matplotlib
         matplotlib.use('Agg')
@@ -266,5 +274,13 @@ def main(args):
     phonopy_qha.write_gruneisen_temperature()
 
 
+def main():
+    if _phonopy_import_error is not None:
+        raise SystemExit(
+            "phonopy-qha requires the optional phonopy dependency: "
+            "pip install 'csld[qha]'\n(%s)" % (_phonopy_import_error,))
+    run(get_options())
+
+
 if __name__ == "__main__":
-    main(get_options())
+    main()
diff --git a/scripts/polaron_main b/csld/cli/polaron_main.py
similarity index 99%
rename from scripts/polaron_main
rename to csld/cli/polaron_main.py
index 37f7d5b..f05ff91 100755
--- a/scripts/polaron_main
+++ b/csld/cli/polaron_main.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python3
+"""``polaron_main`` -- polaron/supercell utility tasks."""
 
 import sys
 import numpy as np
@@ -22,6 +22,11 @@
 
 debug_level = 1
 
+#: Parsed command line options, bound by main(). Read as a global by task1(),
+#: standardize_cell() and cuboid_sites().
+options = None
+
+
 def init_cmdline_settingfile():
     """
     Parse command line options
@@ -971,7 +976,11 @@ def estimate_num_supercell(prim, scmat, cutoff):
     return nmin, nrec 
 
 
-if __name__ == '__main__':
+def main():
+    # task1(), standardize_cell() and cuboid_sites() read `options` as a module
+    # global, so keep it bound at module scope as it was when this file was a
+    # top-level script.
+    global options
     options = init_cmdline_settingfile()
 
     if options.task == '0':
@@ -1359,3 +1368,7 @@ def estimate_num_supercell(prim, scmat, cutoff):
     else:
         print("Unknown task", options.task)
         exit(-1)
+
+
+if __name__ == '__main__':
+    main()
diff --git a/csld/lattice_dynamics.py b/csld/lattice_dynamics.py
index 0c4b93c..7ebea3c 100644
--- a/csld/lattice_dynamics.py
+++ b/csld/lattice_dynamics.py
@@ -19,7 +19,7 @@
 from .structure import SupercellStructure
 from .util.mathtool import MixedIndexForImproper, IntegerDigits, tensor_constraint, RMS, mychop
 from .util.tool import pad_right, matrix2text
-from _c_util import fct_trans_c, ld_get_correlation, get_nullspace, init_ldff_basis, ldff_get_corr
+from csld._c_util import fct_trans_c, ld_get_correlation, get_nullspace, init_ldff_basis, ldff_get_corr
 from .coord_utils import ReadPBC2Cart
 from .util.string_utils import str2arr, str2bool
 
@@ -583,7 +583,7 @@ def save_fcshengbte_original(self, sol, ord, tol=1e-20, output_ijkl=True):
         assert ord in (3,4), "Only order 3 or 4 FCTs are accepted by shengbte, got %d"%(ord)
         import io, re
         import os.path
-        from f_util import f_util
+        from csld.f_util import f_util
 #        uscale = u_list[0]
         np.savetxt('sol',sol)
 #        print('SOL : ',len(sol),'\n',sol)
@@ -770,7 +770,7 @@ def get_hessian_dipole(s):
         """
         :param s:  structure with epsilon_inf and born_charge
         """
-        from f_phonon import f_phonon
+        from csld.phonon.f_phonon import f_phonon
         if s.intensive_properties['epsilon_inf'] is None:
             return np.zeros((3*s.num_sites,3*s.num_sites))
         return f_phonon.get_fcm_dipole(s.lattice.matrix.T, s.lattice.inv_matrix, 1E-18, s.cart_coords.T,
@@ -781,7 +781,7 @@ def get_hessian_dipole_corrected(self, s):
         """
         s: supercell
         """
-        from f_phonon import f_phonon
+        from csld.phonon.f_phonon import f_phonon
         fcm_dp = self.get_hessian_dipole(s)
         if self._dpcor is None:
             return fcm_dp
@@ -826,7 +826,7 @@ def pairinfo_to_supercell(prim, sc, ijk_prim, typ_prim, fcm_prim):
 
 
     def get_dpcor(self, bondlen, errtol=1e-7):
-        from f_phonon import f_phonon
+        from csld.phonon.f_phonon import f_phonon
         offd=[[0,0,1],[1,2,2]]
         offdflatUp = [1,2,5]
         offdflatDn = [3,6,7]
diff --git a/csld/meson.build b/csld/meson.build
new file mode 100644
index 0000000..bb12682
--- /dev/null
+++ b/csld/meson.build
@@ -0,0 +1,28 @@
+py.install_sources(
+  [
+    '__init__.py',
+    'analyzer.py',
+    'atomic_model.py',
+    'basic_lattice_model.py',
+    'cluster.py',
+    'common_main.py',
+    'composition.py',
+    'config.py',
+    'coord_utils.py',
+    'csld_main_functions.py',
+    'interface_vasp.py',
+    'lattice.py',
+    'lattice_dynamics.py',
+    'operations.py',
+    'phonon_mode.py',
+    'sites.py',
+    'structure.py',
+    'symm_kpts.py',
+    'symmetry_structure.py',
+  ],
+  subdir: 'csld',
+)
+
+subdir('cli')
+subdir('phonon')
+subdir('util')
diff --git a/csld/phonon/meson.build b/csld/phonon/meson.build
new file mode 100644
index 0000000..3f4c0b7
--- /dev/null
+++ b/csld/phonon/meson.build
@@ -0,0 +1,49 @@
+# f2py extension, installed as csld.phonon.f_phonon
+#
+# The free-form dostet.f90 / pdstet.f90 / setk06.f90 supersede the fixed-form
+# .f files of the same name, which are dead and must not be compiled as well.
+f_phonon_sources = [
+  'f_phonon.f90',
+  # tetrahedron method
+  'dostet.f90',
+  'pdstet.f90',
+  'setk06.f90',
+  # Ewald / dipole-dipole (fixed form)
+  'DM_dipole_dipole.f',
+  'gbox.f',
+  'tripl.f',
+  'derfc.f',
+  'find_Ewald_eta_screened.f',
+  'DMstandard.f',
+  'lattc.f',
+  'cross.f',
+  'Fewald.f',
+]
+
+f_phonon_gen = custom_target(
+  'f_phononmodule.c',
+  input: ['f_phonon.f90'],
+  output: ['f_phononmodule.c', 'f_phonon-f2pywrappers2.f90'],
+  command: [f2py, '@INPUT@', '-m', 'f_phonon', '--lower', '--build-dir', '@OUTDIR@'],
+)
+
+py.extension_module(
+  'f_phonon',
+  f_phonon_sources + [f_phonon_gen, fortranobject_c],
+  include_directories: [inc_np, inc_f2py],
+  dependencies: [py_dep, lapack],
+  fortran_args: fortran_args,
+  install: true,
+  subdir: 'csld/phonon',
+)
+
+py.install_sources(
+  [
+    '__init__.py',
+    'bandstructure.py',
+    'dos.py',
+    'phonon.py',
+    'plotter.py',
+  ],
+  subdir: 'csld/phonon',
+)
diff --git a/csld/phonon/phonon.py b/csld/phonon/phonon.py
index b77e77e..a2638a2 100644
--- a/csld/phonon/phonon.py
+++ b/csld/phonon/phonon.py
@@ -27,7 +27,7 @@
     plt = None
     pass
 
-from f_phonon import f_phonon
+from csld.phonon.f_phonon import f_phonon
 
 #import logging
 #logger = logging.getLogger(__name__)
diff --git a/csld/structure.py b/csld/structure.py
index 4cffdf4..1e26d2b 100644
--- a/csld/structure.py
+++ b/csld/structure.py
@@ -1739,7 +1739,7 @@ def set_coords(self, c, cart=True):
             self[i].set_coords(c[i], cart=cart)
 
     def get_order_wrt(self, p1, inverse=False, tol=1E-4):
-        from _c_util import get_structure_ordering
+        from csld._c_util import get_structure_ordering
         if p1 is None:
             return list(range(self.num_sites))
 
@@ -1803,7 +1803,7 @@ def match(self, p2, tol_match=0.15, tol_distinct=1.0):
         self: ideal (supercell) structure
         p2: structure with distortion, defect and/or disorder
         """
-        from f_util import f_util
+        from csld.f_util import f_util
         return f_util.match_structure(self.lattice.matrix.T, self.frac_coords.T, self.atomic_numbers, p2.frac_coords.T, p2.atomic_numbers, tol_match, tol_distinct)
 
 
diff --git a/csld/symmetry_structure.py b/csld/symmetry_structure.py
index 94809b7..efb8921 100644
--- a/csld/symmetry_structure.py
+++ b/csld/symmetry_structure.py
@@ -12,7 +12,7 @@
 from .structure import Structure
 try:
     from .util.mathtool import tensor_constraint, get_symmetrized_lsq, mychop
-    from _c_util import fct_trans_c
+    from csld._c_util import fct_trans_c
 except:
     pass
 
@@ -52,7 +52,7 @@ def __init__(self, structure, spacegroup, equivalent_positions, syminfo, tol=1E-
         self.orbit_of_l = inv
         self.reprL_of_l = equivalent_positions
         self._l_list_of_orb = None # wait till setup of orbits
-        self.wyckoffs = syminfo.get_symmetry_dataset()['wyckoffs']
+        self.wyckoffs = syminfo.get_symmetry_dataset().wyckoffs
         # print("DEBUG equiv", equivalent_positions,'spg=',spacegroup,spacegroup.__class__)
         # print("debug u=", u)
         # print("debug inv=", inv)
diff --git a/csld/util/io_utils.py b/csld/util/io_utils.py
index a0b66ef..7a53171 100644
--- a/csld/util/io_utils.py
+++ b/csld/util/io_utils.py
@@ -193,7 +193,7 @@ def clean_json(input_json, strict=False):
         if not strict:
             return str(input_json)
         else:
-            if isinstance(input_json, basestring):
+            if isinstance(input_json, str):
                 return str(input_json)
             else:
                 return clean_json(input_json.to_dict, strict=strict)
diff --git a/csld/util/mathtool.py b/csld/util/mathtool.py
index 48556e5..5d9b3aa 100644
--- a/csld/util/mathtool.py
+++ b/csld/util/mathtool.py
@@ -5,7 +5,7 @@
 from scipy.sparse import lil_matrix as spmat
 # from .rref import rref
 #from .cy_rref import cy_rref
-from _c_util import fct_trans_c, get_nullspace
+from csld._c_util import fct_trans_c, get_nullspace
 
 
 def cofactor(M):
diff --git a/csld/util/meson.build b/csld/util/meson.build
new file mode 100644
index 0000000..9a1c3fb
--- /dev/null
+++ b/csld/util/meson.build
@@ -0,0 +1,23 @@
+py.install_sources(
+  [
+    '__init__.py',
+    'bandstructure.py',
+    'io_utils.py',
+    'mathtool.py',
+    'periodic_table.py',
+    'physical_constants.py',
+    'rref.py',
+    'string_utils.py',
+    'tool.py',
+    'tool_for_original_shengbte.py',
+    'units.py',
+  ],
+  subdir: 'csld/util',
+)
+
+# periodic_table.py loads this with os.path.dirname(__file__), so it has to sit
+# next to the module.
+install_data(
+  'periodic_table.json',
+  install_dir: py.get_install_dir() / 'csld/util',
+)
diff --git a/cssolve/bregman_func.py b/cssolve/bregman_func.py
index cce1c96..5cc4697 100755
--- a/cssolve/bregman_func.py
+++ b/cssolve/bregman_func.py
@@ -15,7 +15,7 @@
 except:
     from scipy.sparse.linalg.eigen.arpack import eigsh
 try:
-    from bregman import bregman
+    from cssolve.bregman import bregman
 except ImportError:
     print("Failed to import Bregman subroutines")
     pass
diff --git a/cssolve/csfit.py b/cssolve/csfit.py
index f026ec5..ef382c7 100755
--- a/cssolve/csfit.py
+++ b/cssolve/csfit.py
@@ -16,7 +16,7 @@
     print("Failed to import bregman_func")
     pass
 try:
-    from bcs_driver import bcs_driver
+    from cssolve.bcs_driver import bcs_driver
 except ImportError:
     #print("Failed to import bcs solver")
     pass
diff --git a/cssolve/meson.build b/cssolve/meson.build
new file mode 100644
index 0000000..caf03a3
--- /dev/null
+++ b/cssolve/meson.build
@@ -0,0 +1,26 @@
+# f2py extension, installed as cssolve.bregman
+bregman_gen = custom_target(
+  'bregmanmodule.c',
+  input: ['bregman.f90'],
+  output: ['bregmanmodule.c', 'bregman-f2pywrappers2.f90'],
+  command: [f2py, '@INPUT@', '-m', 'bregman', '--lower', '--build-dir', '@OUTDIR@'],
+)
+
+py.extension_module(
+  'bregman',
+  ['bregman.f90', bregman_gen, fortranobject_c],
+  include_directories: [inc_np, inc_f2py],
+  dependencies: [py_dep, lapack],
+  fortran_args: fortran_args,
+  install: true,
+  subdir: 'cssolve',
+)
+
+py.install_sources(
+  [
+    '__init__.py',
+    'bregman_func.py',
+    'csfit.py',
+  ],
+  subdir: 'cssolve',
+)
diff --git a/cssolve/setup.py b/cssolve/setup.py
deleted file mode 100644
index 38616c8..0000000
--- a/cssolve/setup.py
+++ /dev/null
@@ -1,31 +0,0 @@
-#!/usr/bin/env python3
-
-import glob
-import os
-import subprocess
-
-#from ez_setup import use_setuptools
-#use_setuptools()
-from setuptools import setup, find_packages, Extension
-from Cython.Distutils import build_ext
-from Cython.Build import cythonize
-import numpy
-
-
-subprocess.call(['make'])
-setup(
-     name="cssolve",
-     packages=find_packages("cssolve"),
-     version="0.7",
-     install_requires=["numpy>=1.9", "scipy>=0.13", "matplotlib>=1.4"],
-     scripts="csfit.py",
-     package_data={'cssolve': ['bregman*.so']},
-#     license="MIT",
-     description="CS solver for CSLD package"
-)
-
-#from numpy.distutils.core import setup, Extension
-#setup(ext_modules= [Extension('bregman', 
-#    sources=['csld/cs_fitting/bregman.f90'],
-#    extra_f90_compile_args=["-cpp", "-heap-arrays"])])
-
diff --git a/environment.yml b/environment.yml
new file mode 100644
index 0000000..7846210
--- /dev/null
+++ b/environment.yml
@@ -0,0 +1,41 @@
+# Conda environment for building and running CSLD.
+#
+#   conda env create -f environment.yml
+#   conda activate csld
+#   pip install .
+#
+# The `compilers` metapackage brings in matching C, C++ and Fortran compilers,
+# so no system toolchain is required. Activating the environment is what puts
+# meson/ninja on PATH and sets CC/CXX/FC, so do that before `pip install`.
+name: csld
+channels:
+  - conda-forge
+dependencies:
+  # Pin a conventional (GIL-enabled) interpreter. CSLD also builds and runs on
+  # newer versions, including the free-threaded 3.14 build, but the extensions
+  # do not declare free-threading support so the GIL is re-enabled at import.
+  - python=3.12
+
+  # toolchain
+  - compilers          # c-compiler, cxx-compiler, fortran-compiler
+  - libblas
+  - liblapack          # no .pc file on conda-forge; meson falls back to -llapack
+  - libopenblas
+  - pkg-config         # used by meson's dependency() lookups
+
+  # build backend
+  - meson>=1.2
+  - meson-python>=0.15
+  - ninja
+  - cython>=3.0
+
+  # runtime
+  - numpy>=1.23
+  - scipy>=1.9
+  - matplotlib-base>=3.5
+  - spglib>=2.5
+
+  # optional / development
+  - scikit-learn
+  - pytest>=7
+  - pip
diff --git a/examples/Si/fit-forcefield.sh b/examples/Si/fit-forcefield.sh
index 66e190e..6a94707 100644
--- a/examples/Si/fit-forcefield.sh
+++ b/examples/Si/fit-forcefield.sh
@@ -1,4 +1,4 @@
-SCRIPT=../../scripts/csld_main
+SCRIPT=csld_main
 # first fit force field coefficients only
 $SCRIPT --ldff_step 2 -f csld.in-forcefield --phonon_step -1 
 # next read coefficients and fit FCT only
diff --git a/manual.rst b/manual.rst
index 7807abf..defbca5 100644
--- a/manual.rst
+++ b/manual.rst
@@ -31,48 +31,64 @@ Installation
 
 Hereby we refer to all files from the main source code directory. User input commands start with "$", followed by output without leading "$"::
 
- $ ls setup.py
- setup.py
+ $ ls pyproject.toml
+ pyproject.toml
 
 - Prerequisite
 
-  - Python 3
-  - Python packages: numpy, scipy, matplotlib, ConfigParser
-  - Library `spglib`_
-  - f2py3 associated with Python3. If it is available in another name, e.g. f2py, make it available by e.g. alias f2py3=f2py
+  - Python 3.9 or newer, and ``pip``
   - C++ and Fortran 90 compilers
+  - BLAS/LAPACK
+
+  The Python dependencies (numpy, scipy, matplotlib and `spglib`_) are installed automatically by ``pip``; so are the build tools (`meson-python`_, meson, ninja and Cython).
 
   If the prerequisites are not met, please either ask your sys admin to install them, or install your own python 3 environment, e.g. `miniconda`_ or `anaconda`_. The codes were tested on conda installations on Linux and MacOS.
 
-.. _spglib: https://atztogo.github.io/spglib/
+.. _spglib: https://spglib.readthedocs.io/
+.. _meson-python: https://mesonbuild.com/meson-python/
 .. _miniconda: https://docs.conda.io/en/latest/miniconda.html
 .. _anaconda: https://docs.anaconda.com/anaconda/install/
 
-- If you have no plan to modify/contribute to the source codes, go to the code directory and install for all users
+- To install, go to the code directory and run
 
 ::
 
-   $ python3 setup.py install
+   $ pip install .
 
 - Install for yourself without admin rights
 
 ::
 
-   $ python3 setup.py install --user
+   $ pip install --user .
+
+- If you do not have compilers or BLAS/LAPACK on your system, the bundled conda environment supplies them. Activating the environment is what puts the compilers and the meson build tools on your ``PATH``, so do that before installing.
+
+::
+
+   $ conda env create -f environment.yml
+   $ conda activate csld
+   $ pip install .
+
+  On Debian/Ubuntu the system packages are ``gfortran liblapack-dev libblas-dev``; on MacOS, ``brew install gcc lapack``.
 
-- If you plan to develop (i.e. modify) the codes
+- If you plan to develop (i.e. modify) the codes, install in editable mode. The Python modules are then served straight from your checkout, and the compiled extensions are rebuilt automatically when they change.
 
 ::
 
-  $ python3 setup.py develop
+  $ pip install -e . --no-build-isolation
 
-- Develop without admin rights
+- To verify the installation
 
 ::
 
-  $ python3 setup.py develop --user
+  $ pip install '.[test]'
+  $ pytest
+
+  This runs the ``examples/Si`` and ``examples/NaCl`` calculations end to end and checks the resulting phonon spectra.
+
+- After installation the commands ``csld_main``, ``cs-fit``, ``polaron_main``, ``phonopy-qha``, ``get-force.sh`` and ``prepare-volume.sh`` are available on your ``PATH``. ``phonopy-qha`` additionally needs `phonopy`_, installed with ``pip install '.[qha]'``.
 
-- If compilation fails, you may need to modify the Makefile manually.
+.. _phonopy: https://phonopy.github.io/phonopy/
 
 ***********
 How to cite
diff --git a/meson.build b/meson.build
new file mode 100644
index 0000000..5fa5b4f
--- /dev/null
+++ b/meson.build
@@ -0,0 +1,85 @@
+project(
+  'csld',
+  'c', 'cpp', 'fortran',
+  version: '1.1.0',
+  license: 'MIT',
+  meson_version: '>=1.2.0',
+  default_options: [
+    'buildtype=release',
+    'cpp_std=c++17',
+  ],
+)
+
+add_languages('cython', native: false)
+
+py = import('python').find_installation(pure: false)
+py_dep = py.dependency()
+
+fc = meson.get_compiler('fortran')
+cpp = meson.get_compiler('cpp')
+
+# ---------------------------------------------------------------------------
+# numpy / f2py
+# ---------------------------------------------------------------------------
+incdir_numpy = run_command(
+  py, ['-c', 'import numpy; print(numpy.get_include())'],
+  check: true,
+).stdout().strip()
+
+incdir_f2py = run_command(
+  py, ['-c', 'from numpy import f2py; print(f2py.get_include())'],
+  check: true,
+).stdout().strip()
+
+inc_np = include_directories(incdir_numpy)
+inc_f2py = include_directories(incdir_f2py)
+
+# Every f2py extension must link the f2py runtime shipped inside numpy.
+fortranobject_c = incdir_f2py / 'fortranobject.c'
+
+f2py = [py, '-m', 'numpy.f2py']
+
+# ---------------------------------------------------------------------------
+# BLAS / LAPACK
+#
+# cssolve/bregman.f90 uses dgemm + dsyev; csld/phonon uses ZHEEV and DGEEV.
+# Try pkg-config first (conda-forge and Debian both ship lapack.pc /
+# openblas.pc), then fall back to a plain library search.
+# ---------------------------------------------------------------------------
+lapack = dependency('lapack', required: false)
+if not lapack.found()
+  lapack = dependency('openblas', required: false)
+endif
+if not lapack.found()
+  lapack = fc.find_library('lapack', required: false)
+endif
+if not lapack.found()
+  lapack = fc.find_library('openblas', required: true)
+endif
+
+# ---------------------------------------------------------------------------
+# Fortran flags
+#
+# -cpp: bregman.f90, f_phonon.f90 and several of the fixed-form .f files
+#       contain preprocessor directives (the old Makefile passed -cpp globally).
+# -fallow-argument-mismatch: gfortran >= 10 rejects the rank/type mismatches in
+#       the legacy fixed-form sources by default.
+# ---------------------------------------------------------------------------
+_fflags = ['-cpp']
+fortran_args = _fflags + fc.get_supported_arguments(
+  '-fallow-argument-mismatch',
+  '-ffree-line-length-none',
+)
+
+subdir('compile/c_util')
+subdir('compile/f_util')
+subdir('cssolve')
+subdir('csld')
+
+# Shell helpers, installed next to the console scripts.
+install_data(
+  'scripts/get-force.sh',
+  'scripts/prepare-volume.sh',
+  install_dir: get_option('bindir'),
+  install_mode: 'rwxr-xr-x',
+)
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..2d61a96
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,57 @@
+[build-system]
+build-backend = "mesonpy"
+requires = [
+    "meson-python>=0.15",
+    "meson>=1.2",
+    "ninja",
+    "Cython>=3.0",
+    "numpy>=2.0",
+]
+
+[project]
+name = "csld"
+dynamic = ["version"]
+description = "Compressive sensing lattice dynamics"
+readme = "README.md"
+license = { text = "MIT" }
+authors = [{ name = "Fei Zhou", email = "fei.fzhou@gmail.com" }]
+requires-python = ">=3.9"
+keywords = ["lattice dynamics", "phonons", "compressive sensing", "force constants"]
+classifiers = [
+    "Development Status :: 5 - Production/Stable",
+    "Intended Audience :: Science/Research",
+    "License :: OSI Approved :: MIT License",
+    "Programming Language :: Python :: 3",
+    "Programming Language :: Fortran",
+    "Programming Language :: C++",
+    "Topic :: Scientific/Engineering :: Physics",
+]
+dependencies = [
+    "numpy>=1.23",
+    "scipy>=1.9",
+    "matplotlib>=3.5",
+    "spglib>=2.5",
+]
+
+[project.optional-dependencies]
+# scripts/phonopy-qha -> csld.cli.phonopy_qha
+qha = ["phonopy>=2.7"]
+# optional solver branch in cssolve/csfit.py
+fit = ["scikit-learn"]
+test = ["pytest>=7"]
+
+[project.scripts]
+csld_main = "csld.cli.csld_main:main"
+cs-fit = "csld.cli.cs_fit:main"
+polaron_main = "csld.cli.polaron_main:main"
+phonopy-qha = "csld.cli.phonopy_qha:main"
+
+[project.urls]
+Homepage = "https://github.com/LLNL/csld"
+Source = "https://github.com/LLNL/csld"
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
+markers = [
+    "slow: end-to-end runs of the examples (deselect with '-m \"not slow\"')",
+]
diff --git a/setup.cfg b/setup.cfg
deleted file mode 100644
index fa0d020..0000000
--- a/setup.cfg
+++ /dev/null
@@ -1,8 +0,0 @@
-[build]
-force = 0
-
-[egg_info]
-tag_build = 
-tag_svn_revision = 0
-tag_date = 0
-
diff --git a/setup.py b/setup.py
deleted file mode 100755
index d1c9fe0..0000000
--- a/setup.py
+++ /dev/null
@@ -1,69 +0,0 @@
-#!/usr/bin/env python3
-
-import glob
-import os
-import subprocess
-
-#from ez_setup import use_setuptools
-#use_setuptools()
-from setuptools import setup, find_packages, Extension
-from Cython.Distutils import build_ext
-from Cython.Build import cythonize
-import numpy
-
-#try:
-#    from numpy.distutils.misc_util import get_numpy_include_dirs
-#except ImportError:
-#    print("numpy.distutils.misc_util cannot be imported. Attempting to "
-#          "install...")
-#    subprocess.call(["easy_install", "numpy"])
-#    from numpy.distutils.misc_util import get_numpy_include_dirs
-
-subprocess.call(['make'])
-
-#def package_files(directory):
-#    paths = []
-#    for (path, directories, filenames) in os.walk(directory):
-#        for filename in filenames:
-#            paths.append(os.path.join('..', path, filename))
-#    return paths
-#data_files = package_files('examples')
-
-setup(
-     name="csld",
-     author="Fei Zhou",
-     author_email="fei.fzhou@gmail.com",
-     url="https://to-be-determined",
-     platforms=['any'],
-     packages=find_packages(),
-     version="1.0",
-     install_requires=["numpy>=1.9", "scipy>=0.13", "matplotlib>=1.4", "spglib>=1.9"],
-     package_data={"csld.util": ["*.json"],
-#       'csld': ['../bregman*.so', '../f_phonon*.so', '../bcs_driver*.so','../f_util*.so','../Makefile', '../css*/*.f90','../compile/f_util/*.f90', '../csld/*/*.f90', '../csld/*/*.f']},
-       'csld': ['../bregman*.so', '../f_phonon*.so', '../f_util*.so','../Makefile', '../css*/*.f90','../compile/f_util/*.f90', '../csld/*/*.f90', '../csld/*/*.f']},
-     #data_files=data_files,
-     license="MIT",
-     description="CSLD",
-     long_description="Compressive sensing lattice dynamics",
-     cmdclass={'build_ext': build_ext},
-     ext_modules=[Extension("_c_util", sources=glob.glob('compile/c_util/*.pyx')+glob.glob('compile/c_util/[a-zA-Z]*.cpp'),
-                            extra_compile_args=['-DNDEBUG','-O3','-g0'], extra_link_args=[],
-                            include_dirs=[numpy.get_include(), 'compile/c_util'], language='c++', libraries = ['stdc++'])],
-     scripts=glob.glob("scripts/*")
-)
-
-# from numpy.distutils.core import setup, Extension
-# print(dir(Extension))
-# setup(ext_modules= [
-#     Extension('bregman', sources=['cssolve/bregman.f90'],
-#               extra_link_args=["-llapack"],
-#               #f2py_options=['--link-lapack_opt'],
-#               extra_f90_compile_args=["-cpp", "-heap-arrays"]),
-#     Extension('f_phonon', sources=['csld/phonon/f_phonon.f90'],
-#               #f2py_options=['--link-lapack_opt'],
-#               extra_f90_compile_args=["-cpp", "-heap-arrays"]),
-#     Extension('bcs_driver', sources=['cssolve/bcs_driver.f90',
-#             'cssolve/num_types.f90', 'cssolve/matrix_sets.f90', 'cssolve/laplace.f90', 'cssolve/bcs.f90'],
-#               extra_f90_compile_args=["-cpp", "-heap-arrays"])
-#     ])
-
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 0000000..342f16a
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,44 @@
+"""Make sure the tests exercise the *installed* csld, not the source tree.
+
+``python -m pytest`` prepends the working directory to ``sys.path``, which would
+let the bare ``csld/`` source directory shadow the installed package. That
+directory has no compiled extensions in it (meson builds them out of tree), so
+every ``csld._c_util`` / ``csld.f_util`` import would fail confusingly.
+"""
+
+import os
+import shutil
+import sys
+
+REPO_ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+EXAMPLES = os.path.join(REPO_ROOT, "examples")
+
+sys.path[:] = [
+    p for p in sys.path
+    if os.path.abspath(p or os.getcwd()) != REPO_ROOT
+]
+
+# An editable install legitimately serves the pure-Python modules from the
+# source tree and the extensions from the meson build directory, so check for
+# the extensions rather than for csld.__file__.
+try:
+    import csld._c_util  # noqa: F401
+except ImportError as exc:  # pragma: no cover - configuration error
+    raise RuntimeError(
+        "csld's compiled extensions are not importable (%s).\n"
+        "Install the package first:  pip install .\n"
+        "or, for a development checkout:  pip install -e . --no-build-isolation"
+        % (exc,)
+    ) from exc
+
+import pytest  # noqa: E402
+
+
+@pytest.fixture
+def example_dir(tmp_path):
+    """Copy one of the bundled examples into a scratch directory."""
+    def _copy(name):
+        dest = tmp_path / name
+        shutil.copytree(os.path.join(EXAMPLES, name), dest)
+        return dest
+    return _copy
diff --git a/tests/test_examples.py b/tests/test_examples.py
new file mode 100644
index 0000000..ef7fae9
--- /dev/null
+++ b/tests/test_examples.py
@@ -0,0 +1,86 @@
+"""End-to-end runs of the bundled examples.
+
+These drive the whole pipeline -- symmetry analysis (spglib), cluster
+enumeration and symmetrisation (``_c_util``, ``f_util``), compressive-sensing
+fitting (``bregman``) and phonon calculation (``f_phonon``) -- through the
+``csld_main`` console script, exactly as a user would.
+
+The fit draws random training/holdout subsets (``np.random.choice`` in
+``cssolve.csfit`` and ``random_seed()`` in ``cssolve/bregman.f90``) and neither
+RNG is seeded, so the fitted force constants differ slightly between runs. The
+assertions below therefore check run-to-run stable physics -- branch count,
+acoustic modes vanishing at Gamma, the known optical frequency, LO-TO splitting
+-- rather than exact numbers.
+"""
+
+import subprocess
+import sys
+
+import numpy as np
+import pytest
+
+pytestmark = pytest.mark.slow
+
+
+def run_csld_main(workdir, *args):
+    """Invoke the csld_main entry point in `workdir` and return its output."""
+    proc = subprocess.run(
+        [sys.executable, "-m", "csld.cli.csld_main", *args],
+        cwd=workdir, capture_output=True, text=True,
+    )
+    assert proc.returncode == 0, (
+        "csld_main failed (exit %d)\n--- stdout ---\n%s\n--- stderr ---\n%s"
+        % (proc.returncode, proc.stdout[-4000:], proc.stderr[-4000:])
+    )
+    return proc.stdout
+
+
+def read_dispersion(workdir):
+    d = np.loadtxt(workdir / "phonon-dispersion.out")
+    assert np.isfinite(d).all(), "non-finite phonon frequencies"
+    # column 0 is the path coordinate, the rest are branches
+    return d[:, 0], d[:, 1:]
+
+
+def test_si_example(example_dir):
+    """Si: 2 atoms/cell -> 6 branches, optical mode at Gamma near 15.5 THz."""
+    work = example_dir("Si")
+    out = run_csld_main(work, "-f", "csld.in")
+
+    assert "Phonon done" in out
+    _, freqs = read_dispersion(work)
+    assert freqs.shape[1] == 6
+
+    # acoustic branches go to zero at Gamma (first point of the path)
+    assert np.abs(freqs[0, :3]).max() < 0.05
+    # highest optical frequency of Si is ~15.5 THz
+    assert 14.5 < freqs.max() < 16.5
+    # a 3x3x3 supercell fit leaves small imaginary acoustic modes near Gamma;
+    # anything large means the force constants are broken
+    assert freqs.min() > -1.0
+
+    for name in ("phonon-total-dos.out", "phonon-partial-dos.out",
+                 "thermal_out.txt", "sol_2nd"):
+        assert (work / name).is_file(), name
+
+
+def test_nacl_example(example_dir):
+    """NaCl: polar, so the Born charges must produce an LO-TO splitting.
+
+    This is the test that exercises the fixed-form Ewald/dipole-dipole Fortran
+    (DM_dipole_dipole.f, find_Ewald_eta_screened.f, ...) linked into f_phonon.
+    """
+    work = example_dir("NaCl")
+    out = run_csld_main(work, "-f", "csld.in")
+
+    assert "Phonon done" in out
+    _, freqs = read_dispersion(work)
+    assert freqs.shape[1] == 6
+
+    assert np.abs(freqs[0, :3]).max() < 0.05
+    assert freqs.min() > -1.0
+
+    # at Gamma the two TO modes are degenerate and the LO mode sits well above
+    to1, to2, lo = freqs[0, 3], freqs[0, 4], freqs[0, 5]
+    assert to1 == pytest.approx(to2, abs=0.05), "TO modes should be degenerate"
+    assert lo > to2 + 1.0, "expected LO-TO splitting from the Born charges"
diff --git a/tests/test_extensions.py b/tests/test_extensions.py
new file mode 100644
index 0000000..8cbbbb6
--- /dev/null
+++ b/tests/test_extensions.py
@@ -0,0 +1,106 @@
+"""Smoke tests for the compiled extension modules.
+
+These prove that the four extensions built, installed to the right place inside
+the ``csld``/``cssolve`` packages, and are callable -- i.e. that the f2py and
+Cython halves of the build are wired up correctly.
+"""
+
+import numpy as np
+import pytest
+import scipy.sparse
+
+
+def test_c_util_importable():
+    import csld._c_util  # noqa: F401
+
+
+def test_f_util_importable():
+    from csld.f_util import f_util  # noqa: F401
+
+
+def test_f_phonon_importable():
+    from csld.phonon.f_phonon import f_phonon  # noqa: F401
+
+
+def test_bregman_importable():
+    from cssolve.bregman import bregman  # noqa: F401
+
+
+def test_fct_trans_c_identity():
+    """The identity rotation must map a rank-2 tensor to the identity matrix."""
+    from csld._c_util import fct_trans_c
+
+    gamma = np.eye(3)
+    gm = fct_trans_c(2, 3, gamma, [0, 1])
+    assert gm.shape == (9, 9)
+    np.testing.assert_allclose(gm.todense(), np.eye(9), atol=1e-12)
+
+
+def test_get_nullspace():
+    """A rank-1 3x3 matrix has a 2-dimensional null space."""
+    from csld._c_util import get_nullspace
+
+    mat = scipy.sparse.csr_matrix(np.array([[1.0, 0.0, 0.0],
+                                            [0.0, 0.0, 0.0],
+                                            [0.0, 0.0, 0.0]]))
+    ns = get_nullspace(mat)
+    assert ns.shape[1] == 3
+    assert ns.shape[0] == 2
+    # every returned vector must actually be annihilated by mat
+    np.testing.assert_allclose(mat.dot(ns.todense().T), 0.0, atol=1e-12)
+
+
+def test_get_structure_ordering_permutation():
+    from csld._c_util import get_structure_ordering
+
+    p1 = np.array([[0.0, 0.0, 0.0], [0.5, 0.5, 0.0], [0.5, 0.0, 0.5]])
+    p2 = p1[[2, 0, 1]].copy()
+    order = get_structure_ordering(p1, p2)
+    np.testing.assert_array_equal(np.sort(order), np.arange(3))
+    np.testing.assert_allclose(p2[order], p1, atol=1e-8)
+
+
+def test_f_util_periodic_distance():
+    """Nearest-image distance in a cubic 4 A cell."""
+    from csld.f_util import f_util
+
+    latt = np.eye(3) * 4.0
+    # 0.9 in fractional coords wraps round to -0.1, i.e. 0.4 Angstrom
+    assert f_util.periodic_distance(latt.T, [0.9, 0.0, 0.0]) == pytest.approx(0.4)
+    assert f_util.periodic_distance(latt.T, [0.5, 0.0, 0.0]) == pytest.approx(2.0)
+
+
+def test_bregman_recovers_sparse_signal():
+    """The L1 solver must recover a 3-sparse vector from an underdetermined A.
+
+    This also proves the LAPACK/BLAS link is live -- BregmanFPC calls dgemm and
+    dsyev via the preconditioner path.
+    """
+    from cssolve.bregman import bregman
+
+    rng = np.random.default_rng(0)
+    n_meas, n_var = 40, 60
+    a = rng.standard_normal((n_meas, n_var))
+    u_true = np.zeros(n_var)
+    u_true[[3, 17, 42]] = [2.0, -1.5, 3.0]
+    b = a.dot(u_true)
+
+    u = bregman.bregmanfpc(10, 1e-6, 1e-3, 1e-4, a, b)
+
+    assert u.shape == (n_var,)
+    assert np.isfinite(u).all()
+    # residual should be small and the support should be recovered
+    assert np.linalg.norm(a.dot(u) - b) / np.linalg.norm(b) < 1e-2
+    assert set(np.argsort(np.abs(u))[-3:]) == {3, 17, 42}
+
+
+def test_f_phonon_exposes_expected_routines():
+    """f_phonon's subroutines carry module state, so only check the interface.
+
+    The numerics are covered end to end by tests/test_examples.py.
+    """
+    from csld.phonon.f_phonon import f_phonon
+
+    for name in ('init', 'get_dispersion', 'get_dm', 'calc_thermal',
+                 'init_nac', 'get_fcm_dipole', 'set_dos_en'):
+        assert callable(getattr(f_phonon, name)), name