diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 7ec3738..4ff3bdb 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -25,8 +25,9 @@ jobs:
           - name: Linux (CUDA 12)
     steps:
     - name: Check out
-      uses: actions/checkout@v2
+      uses: actions/checkout@v4
       with:
+        fetch-depth: 0
         submodules: recursive
         
     - name: Show dependency file
@@ -87,7 +88,8 @@ jobs:
         micromamba activate libmobility
         # install conda-build
         micromamba install -c conda-forge conda-build
-        conda build conda
+        export CUDA_VERSION=12
+        conda build devtools/conda-build
       shell: bash -el {0}
 
     - name: Install via pip
diff --git a/.github/workflows/conda.yaml b/.github/workflows/conda.yaml
new file mode 100644
index 0000000..f2263ef
--- /dev/null
+++ b/.github/workflows/conda.yaml
@@ -0,0 +1,44 @@
+name: Build and upload conda packages
+
+on:
+  release:
+    types: ['released', 'prereleased']
+  # lets add on PR for testing
+  # pull_request:
+  #   types: ['opened', 'edited', 'reopened', 'synchronize']
+
+  workflow_dispatch:        # Un comment line if you also want to trigger action manually
+
+jobs:
+  conda_deployment_with_new_tag:
+    name: Conda deployment of package with Python ${{ matrix.python-version }}
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: ["3.9", "3.10", "3.11", "3.12"]
+        cuda-version: ["12.1", "12.6"]
+    steps:
+      - uses: actions/checkout@v4
+        with:
+          fetch-depth: 0
+          submodules: recursive
+      - name: Conda environment creation and activation
+        uses: conda-incubator/setup-miniconda@v2
+        with:
+          python-version: ${{ matrix.python-version }}
+          environment-file: devtools/conda-envs/build_env.yaml    # Path to the build conda environment
+          auto-update-conda: false
+          auto-activate-base: false
+          show-channel-urls: true
+      - name: Build and upload the conda packages
+        uses: uibcdf/action-build-and-upload-conda-packages@v1.3.0
+        # Export cuda-version as an environment variable
+        env:
+          CUDA_VERSION: ${{ matrix.cuda-version }}
+        with:
+          meta_yaml_dir: devtools/conda-build
+          python-version: ${{ matrix.python-version }} # Values previously defined in `matrix`
+          user: stochasticHydroTools
+          label: auto
+          overwrite: true
+          token: ${{ secrets.ANACONDA_TOKEN }}
diff --git a/README.md b/README.md
index 96d5515..2be978e 100644
--- a/README.md
+++ b/README.md
@@ -25,7 +25,7 @@ This repository is organized into the following directories:
 
 - **include/**: Includes the essential C++ base classes and utility files needed to construct the modules.  
 
-- **conda/**: Contains a meta.yaml file to build a conda package for libMobility using conda-build.  
+- **devtools/**: Contains dev specific scripts and files, like a meta.yaml file to build a conda package for libMobility using conda-build.  
 
 - **docs/**: Contains the source files for the documentation.  
 
@@ -70,12 +70,23 @@ The solvers constructor will check the provided configuration and throw an error
 
 
 ## How to use this repo
+
+### Installation
+
+#### Installing a pre-built binary
+
+You can install libMobility latest release through our conda channel:
+
+```shell
+$ conda install -c conda-forge -c stochasticHydroTools libmobility
+```
+
+#### Compiling from source
+
 Be sure to clone this repository recursively (using ```git clone --recurse```).  
 
 After compilation (see below) you will have all the tools mentioned above available for each solver.
 
-## Compilation
-
 We recommend working with a [conda](https://docs.conda.io/en/latest/) environment. The file environment.yml contains the necessary dependencies to compile and use the library.
 
 You can create the environment with:
@@ -108,7 +119,7 @@ The following variables are available to customize the compilation process:
   * DOUBLEPRECISION : If this variable is defined libMobility is compiled in double precision (single by default).  
 
 
-## Python Usage
+### Python Usage
 
 Importing libMobility will make available any module under "solvers".  
 
@@ -119,12 +130,13 @@ Calling
 ```
 will provide more in depth information about the solver.  
 
-## C++ Usage
+### C++ Usage
 
 In order to use a module called SolverName, the header solvers/SolverName/mobility.h must be included.  
 If the module has been compiled correctly the definitions required for the functions in mobility.h will be available at solvers/SolverName/mobility.so.  
 An example is available in cpp/example.cpp.  
-## Adding a new solver
+
+### Adding a new solver
 
 Solvers must be added following the ```libmobility::Mobility``` directives (see include/MobilityInterface). This C++ base class joins every solver under a common interface. When the C++ interface is prepared, the python bindings can be added automagically by using ```pythonify.h``` (a tool under MobilityInterface).  
 
diff --git a/conda/README.md b/devtools/conda-build/README.md
similarity index 100%
rename from conda/README.md
rename to devtools/conda-build/README.md
diff --git a/conda/build.sh b/devtools/conda-build/build.sh
similarity index 100%
rename from conda/build.sh
rename to devtools/conda-build/build.sh
diff --git a/conda/meta.yaml b/devtools/conda-build/meta.yaml
similarity index 64%
rename from conda/meta.yaml
rename to devtools/conda-build/meta.yaml
index 3c10324..a32be2d 100644
--- a/conda/meta.yaml
+++ b/devtools/conda-build/meta.yaml
@@ -1,25 +1,29 @@
 package:
   name: libmobility
-  version: 0.1.0
+  version: {{ GIT_DESCRIBE_TAG }}
 
 source:
-  path: ../  # Path to your project's source code
+  git_url: ../../
+
+build:
+  number: 0
+  string: cuda{{ CUDA_VERSION }}py{{ CONDA_PY }}h{{ PKG_HASH }}_{{ PKG_BUILDNUM }}
     
 requirements:
 
   build:
     - cmake >=3.22
-    - cuda-version 12.*
+    - cuda-version {{ CUDA_VERSION }}
     - gxx
     - cuda-libraries-dev
     - cuda-nvcc 
     - make
     - mkl-devel
     - pybind11
-    - python 3.11.*
+    - python
   run:
-    - python 3.11.*
-    - cuda-version 12.*
+    - python
+    - cuda-version >={{ CUDA_VERSION }}
     - numpy
     - cuda-libraries
     - mkl
@@ -33,7 +37,7 @@ test:
   source_files:
     - tests
   commands:
-    - pytest -vs -k SelfMobility tests/test_*py
+    - pytest -vs -k "SelfMobility and not fluct" tests/test_*py
     
 about:
   home: https://github.com/stochasticHydroTools/libmobility
diff --git a/devtools/conda-envs/build_env.yaml b/devtools/conda-envs/build_env.yaml
new file mode 100644
index 0000000..297f3a2
--- /dev/null
+++ b/devtools/conda-envs/build_env.yaml
@@ -0,0 +1,7 @@
+channels:    # write here the list of channels to look for your library dependencies
+  - conda-forge
+  - default
+
+dependencies:   # Keep this block with only these two packages
+  - anaconda-client
+  - conda-build
diff --git a/docs/source/installation.rst b/docs/source/installation.rst
index b06ea61..6d115a8 100644
--- a/docs/source/installation.rst
+++ b/docs/source/installation.rst
@@ -1,11 +1,32 @@
 Installation
 ============
 
-We recommend working with a `conda <https://docs.conda.io/en/latest/>`_ environment. The file ``environment.yml`` contains the necessary dependencies to compile and use the library.
+We recommend working with a `conda <https://docs.conda.io/en/latest/>`_ environment.
+
+You can install libMobility's latest release through our conda channel:
+
+.. code-block:: shell
+
+    $ conda install -c conda-forge -c stochasticHydroTools libmobility
+
+
+
+Compilation from source
+~~~~~~~~~~~~~~~~~~~~~~~
+
+If you want to compile the library from source, you can clone the repository with:
+
+.. code-block:: shell
+
+    $ git clone --recursive https://github.com/stochasticHydroTools/libMobility
+
+The ``--recursive`` flag is necessary to clone the submodules of the repository.
 
 Getting dependencies
 --------------------
 
+The file ``environment.yml`` contains the necessary dependencies to compile and use the library.
+
 You can create the environment with:
 
 .. code-block:: shell
@@ -20,9 +41,27 @@ Then, activate the environment with:
 
 .. hint:: At the moment we offer several installation methods: via pip, conda or building from source. We recommend using pip or conda first, resorting to building from source only if the other methods do not work for you.
 
+Building from source
+--------------------
+
+CMake is used for compilation under the hood. After installing the dependencies you can compile and install everything with:
+
+.. code-block:: shell
+
+    $ mkdir build && cd build
+    $ cmake -DCMAKE_INSTALL_PREFIX=$CONDA_PREFIX ..
+    $ make all install
+
+It is advisable to install the library in the conda environment, so that the python bindings are available. The environment variable ``$CONDA_PREFIX`` is set to the root of the conda environment.
+
+After compilation, the python bindings will be available in the conda environment under the name ``libMobility``. See :ref:`usage` for more information.
+
+The following variables are available to customize the compilation process:
+
+- ``DOUBLEPRECISION``: If this variable is defined, libMobility is compiled in double precision (single by default).
 	  
-Installing via pip
-------------------
+Alternative: Installing via pip
+-------------------------------
 
 After installing the dependencies, you can install the library with pip. Go to the root of the repository and run:
 
@@ -30,11 +69,11 @@ After installing the dependencies, you can install the library with pip. Go to t
 
     $ pip install .
     
+   
+Alternative: Building a conda package
+-------------------------------------
 
-Building a conda package
-------------------------
-
-Building the conda package only requires the conda-build package. You can install it with:
+Building the conda package only requires the conda-build package (dependencies will be fetched automatically). You can install it with:
 
 .. code-block:: shell
 
@@ -44,7 +83,7 @@ You can build a conda package with the following command from the root of the re
 
 .. code-block:: shell
 		
-    $ conda build conda
+    $ conda build devtools/conda-build
 
 This will build and test the package, which you can install in any environment with:
 
@@ -55,22 +94,3 @@ This will build and test the package, which you can install in any environment w
 Conda will automatically install all the dependencies needed to run the library.
 
 
-Building from source
---------------------
-
-
-CMake is used for compilation, you can compile and install everything with:
-
-.. code-block:: shell
-
-    $ mkdir build && cd build
-    $ cmake -DCMAKE_INSTALL_PREFIX=$CONDA_PREFIX ..
-    $ make all install
-
-It is advisable to install the library in the conda environment, so that the python bindings are available. The environment variable ``$CONDA_PREFIX`` is set to the root of the conda environment.
-
-After compilation, the python bindings will be available in the conda environment under the name ``libMobility``. See :ref:`usage` for more information.
-
-The following variables are available to customize the compilation process:
-
-- ``DOUBLEPRECISION``: If this variable is defined, libMobility is compiled in double precision (single by default).
diff --git a/setup.py b/setup.py
index da965e7..e0e3341 100644
--- a/setup.py
+++ b/setup.py
@@ -4,6 +4,16 @@
 import os
 import sys
 
+try:
+    version = (
+        subprocess.check_output(["git", "describe", "--abbrev=0", "--tags"])
+        .strip()
+        .decode("utf-8")
+    )
+except:
+    print("Failed to retrieve the current version, defaulting to 0")
+    version = "0"
+
 
 class CMakeExtension(Extension):
     def __init__(self, name, sourcedir=""):
@@ -50,7 +60,7 @@ def build_extension(self, ext):
 
 setup(
     name="libMobility",
-    version="0.1.0",
+    version=version,
     packages=find_packages(),
     ext_modules=[CMakeExtension("libMobility")],
     cmdclass=dict(build_ext=CMakeBuild),
diff --git a/solvers/DPStokes/extra/poly_fits.h b/solvers/DPStokes/extra/poly_fits.h
index 9bfe59a..878fdb1 100644
--- a/solvers/DPStokes/extra/poly_fits.h
+++ b/solvers/DPStokes/extra/poly_fits.h
@@ -1,7 +1,25 @@
+/* Ryker Fish 2024. Helpers to compute values of c^{-1} from [1].
+
+The ES kernel has a function c(beta) that provides a relationship between the hydrodynamic radius and the kernel
+parameters alpha and beta. There is no closed form for c, but we can approximate c and c^{-1} with a polynomial fit.
+
+We have rH = h*w*c(beta), where rH is the hydrodynamic radius, h is the grid spacing, and w is the width of the kernel. 
+Thus, we can write c^{-1}(rH/(h*w)) = beta.
+
+References:
+[1] Computing hydrodynamic interactions in confined doubly periodic geometries in linear time.
+A. Hashemi et al. J. Chem. Phys. 158, 154101 (2023) https://doi.org/10.1063/5.0141371
+ */
+
 #include <vector>
   
 namespace dpstokes_polys{
 
+  /* Evaluates a polynomial at x with coefficients in descending order, so the highest order coefficients
+  are at the start of polyCoeffs.
+  e.g. for a polynomial of order n, this computes 
+      polyCoeffs[n+1] + polyCoeffs[n]*x + polyCoeffs[n-2]*x^2 + ... + polyCoeffs[0]*x^n
+  */
   double polyEval(std::vector<double> polyCoeffs, double x){
 
     int order = polyCoeffs.size() - 1;
@@ -15,6 +33,7 @@ namespace dpstokes_polys{
     return accumulator;
   }
 
+  // Coefficients for the polynomial fit of c^{-1} from [1]
   std::vector<double> cbetam_inv = {
   4131643418.193291,
   -10471683395.26777,
diff --git a/solvers/DPStokes/mobility.h b/solvers/DPStokes/mobility.h
index 39348cb..e065012 100644
--- a/solvers/DPStokes/mobility.h
+++ b/solvers/DPStokes/mobility.h
@@ -1,4 +1,8 @@
 /*Raul P. Pelaez 2022. libMobility interface for UAMMD's DPStokes module
+
+References:
+[1] Computing hydrodynamic interactions in confined doubly periodic geometries in linear time.
+A. Hashemi et al. J. Chem. Phys. 158, 154101 (2023) https://doi.org/10.1063/5.0141371
  */
 #ifndef MOBILITY_SELFMOBILITY_H
 #define MOBILITY_SELFMOBILITY_H
@@ -40,6 +44,9 @@ class DPStokes: public libmobility::Mobility{
 
   void setParametersDPStokes(DPStokesParameters i_dppar){
     this->dppar = i_dppar;
+
+    if(this->dppar.Lx != this->dppar.Ly) throw std::runtime_error("[DPStokes] Only square periodic boxes (Lx = Ly) are currently supported.\n");
+    
     dpstokes = std::make_shared<uammd_dpstokes::DPStokesGlue>();
   }
 
@@ -67,21 +74,13 @@ class DPStokes: public libmobility::Mobility{
     this->dppar.alpha = this->dppar.w*0.5;
     this->dppar.tolerance = 1e-6;
 
-    real N_real = this->dppar.Lx/h; // actual N given L and h
-    int N_up = ceil(N_real);
-    int N_down = floor(N_real);
-    int N;
-    // either N_up or N_down will be a multiple of 2. pick the even one for a more FFT friendly grid.
-    if(N_up % 2 == 0){
-      N = N_up;
-    }else{
-      N = N_down;
-    }
+    int N = floor(this->dppar.Lx/h);
+    N += N % 2;
 
     this->dppar.nx = N;
     this->dppar.ny = N;
 
-    // note: only set up for square boxes
+    // note: this part is only configured for square boxes
     if(this->dppar.allowChangingBoxSize){ // adjust box size to suit h
       this->dppar.Lx = N*h;
       this->dppar.Ly = N*h;
@@ -104,14 +103,9 @@ class DPStokes: public libmobility::Mobility{
     // sets chebyshev node spacing at its coarsest (in the middle) to be h
     real nz_actual = M_PI/(asin(h/H)) + 1;
 
-    // pick nearby N such that 2(Nz-1) is FFT friendly
-    N_up = ceil(nz_actual);
-    N_down = floor(nz_actual);
-    if(N_up % 2 == 1){
-      this->dppar.nz = N_up;
-    } else {
-      this->dppar.nz = N_down;
-    }
+    // pick nearby N such that 2(Nz-1) has two factors of 2 and is FFT friendly
+    this->dppar.nz = floor(nz_actual);
+    this->dppar.nz += (int)ceil(nz_actual) % 2;
 
     dpstokes->initialize(dppar, this->numberParticles);
     Mobility::initialize(ipar);
diff --git a/tests/ref/pair_mobility_bw_ref_noimg.mat b/tests/ref/pair_mobility_bw_ref_noimg.npz
similarity index 99%
rename from tests/ref/pair_mobility_bw_ref_noimg.mat
rename to tests/ref/pair_mobility_bw_ref_noimg.npz
index 97b988b..f7757c1 100644
Binary files a/tests/ref/pair_mobility_bw_ref_noimg.mat and b/tests/ref/pair_mobility_bw_ref_noimg.npz differ
diff --git a/tests/ref/pair_mobility_bw_torque.mat b/tests/ref/pair_mobility_bw_torque.npz
similarity index 98%
rename from tests/ref/pair_mobility_bw_torque.mat
rename to tests/ref/pair_mobility_bw_torque.npz
index 630ed82..896be62 100644
Binary files a/tests/ref/pair_mobility_bw_torque.mat and b/tests/ref/pair_mobility_bw_torque.npz differ
diff --git a/tests/ref/pair_mobility_bw_w4.mat b/tests/ref/pair_mobility_bw_w4.npz
similarity index 96%
rename from tests/ref/pair_mobility_bw_w4.mat
rename to tests/ref/pair_mobility_bw_w4.npz
index b1dde98..e8099fc 100644
Binary files a/tests/ref/pair_mobility_bw_w4.mat and b/tests/ref/pair_mobility_bw_w4.npz differ
diff --git a/tests/ref/pair_mobility_nbody_freespace.mat b/tests/ref/pair_mobility_nbody_freespace.npz
similarity index 81%
rename from tests/ref/pair_mobility_nbody_freespace.mat
rename to tests/ref/pair_mobility_nbody_freespace.npz
index e01c99a..ff716be 100644
Binary files a/tests/ref/pair_mobility_nbody_freespace.mat and b/tests/ref/pair_mobility_nbody_freespace.npz differ
diff --git a/tests/ref/pair_mobility_sc_torque.mat b/tests/ref/pair_mobility_sc_torque.npz
similarity index 98%
rename from tests/ref/pair_mobility_sc_torque.mat
rename to tests/ref/pair_mobility_sc_torque.npz
index 24ecf55..43c2c4e 100644
Binary files a/tests/ref/pair_mobility_sc_torque.mat and b/tests/ref/pair_mobility_sc_torque.npz differ
diff --git a/tests/ref/pair_mobility_sc_w4.mat b/tests/ref/pair_mobility_sc_w4.npz
similarity index 96%
rename from tests/ref/pair_mobility_sc_w4.mat
rename to tests/ref/pair_mobility_sc_w4.npz
index 1797e01..82a8a6a 100644
Binary files a/tests/ref/pair_mobility_sc_w4.mat and b/tests/ref/pair_mobility_sc_w4.npz differ
diff --git a/tests/ref/readme.txt b/tests/ref/readme.txt
index cb80a74..1edcfb6 100644
--- a/tests/ref/readme.txt
+++ b/tests/ref/readme.txt
@@ -9,7 +9,7 @@ sc- slit channel
 -------------- Files --------------
 SELF:
 self_mobility_bw_ref: from the DPStokes repository. Periodized NBody using kernels with wall corrections in double precision.
-self_mobility_bw_ref_noimg.mat: from the DPStokes repository. NBody kernels with wall corrections in double precision, but without using periodized images.
+self_mobility_bw_ref_noimg: from the DPStokes repository. NBody kernels with wall corrections in double precision, but without using periodized images.
 self_mobility_bw_w4, self_mobility_sc_w4: generated from this repository in double precision with default parameters. 
 self_mobility_*_torques: generated from this repository in double precision with default parameters.
 
@@ -17,4 +17,4 @@ PAIR:
 pair_mobility_bw_w4, pair_mobility_sc_w4: generated from this repository in double precision with default parameters. 
 pair_mobility_bw_ref_noimg: from the DPStokes repository. NBody kernels with wall corrections in double precision, but without using periodized images.
 pair_mobility_nbody_freespace: no-wall RPY pair mobility. Generated using mobility matrix numba code from the RigidMultiblob repository.
-pair_mobility_*_torques: generated from this repository in double precision with default parameters.
\ No newline at end of file
+pair_mobility_*_torques: generated from this repository in double precision with default parameters.
diff --git a/tests/ref/self_mobility_bw_ref.mat b/tests/ref/self_mobility_bw_ref.npz
similarity index 95%
rename from tests/ref/self_mobility_bw_ref.mat
rename to tests/ref/self_mobility_bw_ref.npz
index 90d68cf..d3c57bd 100644
Binary files a/tests/ref/self_mobility_bw_ref.mat and b/tests/ref/self_mobility_bw_ref.npz differ
diff --git a/tests/ref/self_mobility_bw_ref_noimg.mat b/tests/ref/self_mobility_bw_ref_noimg.npz
similarity index 92%
rename from tests/ref/self_mobility_bw_ref_noimg.mat
rename to tests/ref/self_mobility_bw_ref_noimg.npz
index b7e4d7f..2a65791 100644
Binary files a/tests/ref/self_mobility_bw_ref_noimg.mat and b/tests/ref/self_mobility_bw_ref_noimg.npz differ
diff --git a/tests/ref/self_mobility_bw_torque.mat b/tests/ref/self_mobility_bw_torque.npz
similarity index 91%
rename from tests/ref/self_mobility_bw_torque.mat
rename to tests/ref/self_mobility_bw_torque.npz
index 5f9a6bf..5f1d155 100644
Binary files a/tests/ref/self_mobility_bw_torque.mat and b/tests/ref/self_mobility_bw_torque.npz differ
diff --git a/tests/ref/self_mobility_bw_w4.mat b/tests/ref/self_mobility_bw_w4.npz
similarity index 83%
rename from tests/ref/self_mobility_bw_w4.mat
rename to tests/ref/self_mobility_bw_w4.npz
index 63744a7..665494c 100644
Binary files a/tests/ref/self_mobility_bw_w4.mat and b/tests/ref/self_mobility_bw_w4.npz differ
diff --git a/tests/ref/self_mobility_sc_torque.mat b/tests/ref/self_mobility_sc_torque.npz
similarity index 91%
rename from tests/ref/self_mobility_sc_torque.mat
rename to tests/ref/self_mobility_sc_torque.npz
index acc4bd3..9d6abff 100644
Binary files a/tests/ref/self_mobility_sc_torque.mat and b/tests/ref/self_mobility_sc_torque.npz differ
diff --git a/tests/ref/self_mobility_sc_w4.mat b/tests/ref/self_mobility_sc_w4.npz
similarity index 84%
rename from tests/ref/self_mobility_sc_w4.mat
rename to tests/ref/self_mobility_sc_w4.npz
index 378db78..c0ea116 100644
Binary files a/tests/ref/self_mobility_sc_w4.mat and b/tests/ref/self_mobility_sc_w4.npz differ
diff --git a/tests/test_fluctuation_dissipation.py b/tests/test_fluctuation_dissipation.py
index efcc207..032b314 100644
--- a/tests/test_fluctuation_dissipation.py
+++ b/tests/test_fluctuation_dissipation.py
@@ -69,8 +69,7 @@ def fluctuation_dissipation_KS(M, fluctuation_method, needsTorques):
     ],
 )
 @pytest.mark.parametrize("hydrodynamicRadius", [0.95, 1.12])
-# @pytest.mark.parametrize("numberParticles", [1,2,10])
-@pytest.mark.parametrize("numberParticles", [10])
+@pytest.mark.parametrize("numberParticles", [1,2,10])
 def test_fluctuation_dissipation_linear_displacements(
     Solver, periodicity, hydrodynamicRadius, numberParticles
 ):
diff --git a/tests/test_initialization.py b/tests/test_initialization.py
index f888c4f..c732603 100644
--- a/tests/test_initialization.py
+++ b/tests/test_initialization.py
@@ -24,7 +24,7 @@ def test_periodic_initialization(Solver, periodicity):
 def test_pse_torques_unsupported():
     hydrodynamicRadius = 1.0
     solver = PSE("periodic", "periodic", "periodic")
-    parameters = sane_parameters[Solver.__name__]
+    parameters = sane_parameters[PSE.__name__]
     solver.setParameters(**parameters)
     numberParticles = 1
     with pytest.raises(RuntimeError):
@@ -41,6 +41,13 @@ def test_invalid_throws(Solver):
     with pytest.raises(RuntimeError):
         solver = Solver("periodicasdas", "periodic", "open")
 
+def test_dpstokes_invalid_box():
+    with pytest.raises(RuntimeError):
+        Lx, Ly = 10, 20
+        solver = DPStokes("periodic", "periodic", "single_wall")
+        params = {"dt": 1, "Lx": Lx, "Ly": Ly, "zmin": 0, "zmax": 19.2}
+        solver.setParameters(**params)
+
 @pytest.mark.parametrize(("NBatch", "NperBatch"),
                          [
                             (1, 12),
diff --git a/tests/test_mobility_matrix.py b/tests/test_mobility_matrix.py
index 2cdffef..d08d40e 100644
--- a/tests/test_mobility_matrix.py
+++ b/tests/test_mobility_matrix.py
@@ -41,6 +41,8 @@ def test_mobility_matrix_linear(
     assert M.shape == (3 * numberParticles, 3 * numberParticles)
     assert M.dtype == precision
     sym = M - M.T
+    atol = 5e-5
+    rtol = 1e-7
     assert np.allclose(
         sym, 0.0, rtol=0, atol=5e-5
     ), f"Mobility matrix is not symmetric within 5e-5, max diff: {np.max(np.abs(sym))}"
@@ -132,7 +134,7 @@ def test_self_mobility_linear_nbody(algorithm):
     forces = np.ones(3, dtype=precision)
     result, _ = solver.Mdot(forces)
     m0 = 1.0 / (6 * np.pi * viscosity * hydrodynamicRadius)
-    assert np.allclose(result, m0 * forces, rtol=0, atol=1e-7)
+    assert np.allclose(result, m0 * forces, rtol=1e-7, atol=1e-7)
 
 @pytest.mark.parametrize("algorithm", ["naive", "block", "fast", "advise"])
 def test_self_mobility_angular_nbody(algorithm):
@@ -211,8 +213,8 @@ def test_pair_mobility_angular_nbody(algorithm):
     parameters = {"algorithm": algorithm}
     solver.setParameters(**parameters)
 
-    ref_file = "./ref/pair_mobility_nbody_freespace.mat"
-    ref = scipy.io.loadmat(ref_file)
+    ref_file = "./ref/pair_mobility_nbody_freespace.npz"
+    ref = np.load(ref_file)
     refM = np.array(ref['M']).astype(precision)
     r_vecs = np.array(ref['r_vecs']).astype(precision)
     a = np.array(ref['a']).astype(precision).flatten()
@@ -233,4 +235,4 @@ def test_pair_mobility_angular_nbody(algorithm):
         solver.setPositions(positions)
 
         M = compute_M(solver, 2, True)
-        assert np.allclose(refM[i], M, atol=1e-6)
\ No newline at end of file
+        assert np.allclose(refM[i], M, atol=1e-6)
diff --git a/tests/test_wall_mobility.py b/tests/test_wall_mobility.py
index 3dcbca5..a18eac2 100644
--- a/tests/test_wall_mobility.py
+++ b/tests/test_wall_mobility.py
@@ -1,12 +1,13 @@
 import pytest
 import numpy as np
-import scipy.interpolate
-import scipy.io
 
 from libMobility import DPStokes, NBody
 from utils import compute_M
 
-self_mobility_params = {
+# NOTE: Some of the following tests will only pass if compiled with double precision.
+# This is because the reference data was generated in double precision.
+
+wall_params = {
     "DPStokes": {"dt": 1, "Lx": 76.8, "Ly": 76.8, "zmin": 0, "zmax": 19.2},
     "NBody": {"algorithm": "advise"},
 }
@@ -14,23 +15,21 @@
 @pytest.mark.parametrize(
     ("Solver", "periodicity", "tol", "start_height", "ref_file"),
     [
-        (DPStokes, ("periodic", "periodic", "single_wall"), 5e-3, 4, "self_mobility_bw_ref.mat"),
-        (DPStokes, ("periodic", "periodic", "single_wall"), 1e-5, 0, "self_mobility_bw_w4.mat"),
-        (DPStokes, ("periodic", "periodic", "two_walls"), 1e-5, 0, "self_mobility_sc_w4.mat"),
-        (NBody, ("open", "open", "single_wall"), 1e-6, 1, "self_mobility_bw_ref_noimg.mat")
+        (DPStokes, ("periodic", "periodic", "single_wall"), 5e-3, 4, "self_mobility_bw_ref.npz"),
+        (DPStokes, ("periodic", "periodic", "single_wall"), 1e-6, 0, "self_mobility_bw_w4.npz"),
+        (DPStokes, ("periodic", "periodic", "two_walls"), 1e-6, 0, "self_mobility_sc_w4.npz"),
+        (NBody, ("open", "open", "single_wall"), 1e-6, 1, "self_mobility_bw_ref_noimg.npz")
     ],
 )
 def test_self_mobility_linear(Solver, periodicity, tol, start_height, ref_file):
-    zmax = 19.2
     xymax = 76.8
-    params = self_mobility_params[Solver.__name__]
-
+    params = wall_params[Solver.__name__]
     needsTorque = False
 
     ref_dir = "./ref/"
-    ref = scipy.io.loadmat(ref_dir + ref_file)
+    ref = np.load(ref_dir + ref_file)
     refM = ref['M']
-    refHeights = ref['heights'][0]
+    refHeights = ref['heights'].flatten()
 
     hydrodynamicRadius = 1.0
     eta = 1/4/np.sqrt(np.pi)
@@ -66,37 +65,29 @@ def test_self_mobility_linear(Solver, periodicity, tol, start_height, ref_file):
         M /= normMat
         allM[i] = M
 
-    # uncomment to save datafile for test plots
-    # scipy.io.savemat('./temp/test_data/test_' + ref_file, {'M': allM, 'heights': refHeights})
-
     diags = [np.diag(matrix) for matrix in allM]
     ref_diags = [np.diag(matrix)[0:3] for matrix in refM] # only take diagonal elements from forces
 
-    diff = np.abs([(diag - ref_diag) for diag, ref_diag in zip(diags, ref_diags)])
-
-    avgErr = np.mean(diff)
-
-    assert avgErr < tol, f"Self mobility does not match reference. Average error: {avgErr}"
+    for diag, ref_diag in zip(diags, ref_diags):
+        assert np.allclose(diag, ref_diag, rtol=tol, atol=tol), f"Self mobility does not match reference"
 
 @pytest.mark.parametrize(
     ("Solver", "periodicity", "tol", "ref_file"),
     [
-        (DPStokes, ("periodic", "periodic", "single_wall"), 5e-5, "pair_mobility_bw_w4.mat"),
-        (DPStokes, ("periodic", "periodic", "two_walls"), 5e-5, "pair_mobility_sc_w4.mat"),
-        (NBody, ("open", "open", "single_wall"), 1e-4, "pair_mobility_bw_ref_noimg.mat"),
+        (DPStokes, ("periodic", "periodic", "single_wall"), 1e-6, "pair_mobility_bw_w4.npz"),
+        (DPStokes, ("periodic", "periodic", "two_walls"), 1e-6, "pair_mobility_sc_w4.npz"),
+        (NBody, ("open", "open", "single_wall"), 1e-4, "pair_mobility_bw_ref_noimg.npz"),
     ],
 )
 def test_pair_mobility_linear(Solver, periodicity, ref_file, tol):
-    zmax = 19.2
     xymax = 76.8
-    params = self_mobility_params[Solver.__name__]
-
+    params = wall_params[Solver.__name__]
     needsTorque = False
 
     ref_dir = "./ref/"
-    ref = scipy.io.loadmat(ref_dir + ref_file)
+    ref = np.load(ref_dir + ref_file)
     refM = ref['M']
-    refHeights = ref['heights'].flatten()
+    refHeights = ref['heights']
     nHeights = len(refHeights)
 
     radH = 1.0 # hydrodynamic radius
@@ -132,9 +123,6 @@ def test_pair_mobility_linear(Solver, periodicity, ref_file, tol):
             M /= normMat
             allM[i][k] = M
 
-    # uncomment to save datafile for test plots
-    # scipy.io.savemat('./temp/test_data/test_' + ref_file, {'M': allM, 'heights': refHeights})
-
     for i in range(0, nSeps):
         for k in range(0, nHeights):
             diff = abs(allM[i,k] - refM[i,k][0:6,0:6])
@@ -143,15 +131,15 @@ def test_pair_mobility_linear(Solver, periodicity, ref_file, tol):
 @pytest.mark.parametrize(
     ("Solver", "periodicity", "ref_file"),
     [
-        (DPStokes, ("periodic", "periodic", "single_wall"), "self_mobility_bw_torque.mat"),
-        (DPStokes, ("periodic", "periodic", "two_walls"), "self_mobility_sc_torque.mat"),
-        (NBody, ("open", "open", "single_wall"), "self_mobility_bw_ref_noimg.mat")
+        (DPStokes, ("periodic", "periodic", "single_wall"), "self_mobility_bw_torque.npz"),
+        (DPStokes, ("periodic", "periodic", "two_walls"), "self_mobility_sc_torque.npz"),
+        (NBody, ("open", "open", "single_wall"), "self_mobility_bw_ref_noimg.npz")
     ],
 )
 def test_self_mobility_angular(Solver, periodicity, ref_file):
     zmax = 19.2
     xymax = 76.8
-    params = self_mobility_params[Solver.__name__]
+    params = wall_params[Solver.__name__]
 
     hydrodynamicRadius = 1.0
     eta = 1/4/np.sqrt(np.pi)
@@ -160,9 +148,9 @@ def test_self_mobility_angular(Solver, periodicity, ref_file):
     tol = 1e-6
 
     ref_dir = "./ref/"
-    ref = scipy.io.loadmat(ref_dir + ref_file)
+    ref = np.load(ref_dir + ref_file)
     refM = ref['M']
-    refHeights = ref['heights'][0]
+    refHeights = ref['heights'].flatten()
 
     precision = np.float32 if Solver.precision == "float" else np.float64
 
@@ -187,6 +175,7 @@ def test_self_mobility_angular(Solver, periodicity, ref_file):
 
     allM = np.zeros((nHeights, 6*numberParticles, 6*numberParticles), dtype=precision)
     for i in range(0,nHeights):
+        # breakpoint()
         positions = np.array([[xymax/2, xymax/2, refHeights[i]]], dtype=precision)
         solver.setPositions(positions)
         
@@ -194,9 +183,6 @@ def test_self_mobility_angular(Solver, periodicity, ref_file):
         M /= normMat
         allM[i] = M
 
-    # uncomment to save datafile for test plots
-    # scipy.io.savemat('./temp/test_data/test_' + ref_file, {'M': allM, 'heights': refHeights})
-
     for i in range(0, nHeights):
         diff = abs(allM[i] - refM[i])
         assert np.all(diff < tol)
@@ -204,15 +190,15 @@ def test_self_mobility_angular(Solver, periodicity, ref_file):
 @pytest.mark.parametrize(
     ("Solver", "periodicity", "ref_file"),
     [
-        (DPStokes, ("periodic", "periodic", "single_wall"), "pair_mobility_bw_torque.mat"),
-        (DPStokes, ("periodic", "periodic", "two_walls"), "pair_mobility_sc_torque.mat"),
-        (NBody, ("open", "open", "single_wall"), "pair_mobility_bw_ref_noimg.mat")
+        (DPStokes, ("periodic", "periodic", "single_wall"), "pair_mobility_bw_torque.npz"),
+        (DPStokes, ("periodic", "periodic", "two_walls"), "pair_mobility_sc_torque.npz"),
+        (NBody, ("open", "open", "single_wall"), "pair_mobility_bw_ref_noimg.npz")
     ],
 )
 def test_pair_mobility_angular(Solver, periodicity, ref_file):
     zmax = 19.2
     xymax = 76.8
-    params = self_mobility_params[Solver.__name__]
+    params = wall_params[Solver.__name__]
     hydrodynamicRadius = 1.0
     eta = 1/4/np.sqrt(np.pi)
     needsTorque = True
@@ -220,7 +206,7 @@ def test_pair_mobility_angular(Solver, periodicity, ref_file):
     tol = 1e-6
 
     ref_dir = "./ref/"
-    ref = scipy.io.loadmat(ref_dir + ref_file)
+    ref = np.load(ref_dir + ref_file)
     refM = ref['M']
     refHeights = ref['heights'].flatten()
 
@@ -260,33 +246,9 @@ def test_pair_mobility_angular(Solver, periodicity, ref_file):
             M /= normMat
             allM[i,k] = M
 
-    # uncomment to save datafile for test plots
-    scipy.io.savemat('./temp/test_data/test_' + ref_file, {'M': allM, 'heights': refHeights, 'seps': seps})
-
     for i in range(0, nSeps):
         for k in range(0, nHeights):
             diff = abs(allM[i,k] - refM[i,k])
             temp = diff < tol
             x = temp[0:6,6:12]
-            # breakpoint()
-            assert np.all(diff < tol)
-
-def checkPairComponent(indx, indy, allM, refM, nSeps, tol):
-
-    indx -= 1 # shift from matlab to python indexing
-    indy -= 1
-    for i in range(0, nSeps):
-
-        x = allM[i, :, indx, indy]
-        x_ref = refM[i, :, indx, indy]
-
-        for j in range(0, len(x)):
-
-            diff = x[j]-x_ref[j]
-
-            if x[j] > tol and x_ref[j] > tol and x_ref[j] != 0:
-                err = diff/x_ref[j]
-            else:
-                err = diff
-
-        assert err < tol, f"Pair mobility does not match reference for component {indx+1}, {indy+1}"
\ No newline at end of file
+            assert np.all(diff < tol)
\ No newline at end of file
diff --git a/tests/utils.py b/tests/utils.py
index 40bc048..fd62749 100644
--- a/tests/utils.py
+++ b/tests/utils.py
@@ -4,7 +4,7 @@
 sane_parameters = {
     "PSE": {"psi": 1.0, "Lx": 32, "Ly": 32, "Lz": 32, "shearStrain": 0.0},
     "NBody": {"algorithm": "advise"},
-    "DPStokes": {"dt": 1, "Lx": 16, "Ly": 16, "zmin": -6, "zmax": 6},
+    "DPStokes": {"dt": 1, "Lx": 16, "Ly": 16, "zmin": -6, "zmax": 6, "allowChangingBoxSize": False},
     "SelfMobility": {"parameter": 5.0},
 }