diff --git a/Changelog.md b/Changelog.md
--- a/Changelog.md
+++ b/Changelog.md
@@ -1,5 +1,9 @@
 # Changelog
 
+## 0.0.2
+  - switching from `PSQueue` to `psqueues` as the former is not maintained
+  - raise upper bounds of base, allows GHC 9.0
+
 ## 0.0.1
   - raising upper bounds of optics and attoparsec
 
diff --git a/ConClusion.cabal b/ConClusion.cabal
--- a/ConClusion.cabal
+++ b/ConClusion.cabal
@@ -5,7 +5,7 @@
 -- see: https://github.com/sol/hpack
 
 name:           ConClusion
-version:        0.0.1
+version:        0.0.2
 synopsis:       Cluster algorithms, PCA, and chemical conformere analysis
 description:    Please see the README on GitLab at <https://gitlab.com/theoretical-chemistry-jena/quantum-chemistry/ConfoCluster>
 category:       Statistics, Chemistry
@@ -60,14 +60,14 @@
       NamedFieldPuns
   ghc-options: -Wall -Wno-unused-top-binds -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wpartial-fields -Wredundant-constraints
   build-depends:
-      PSQueue >=1.1.0.1 && <1.2
-    , aeson ==1.5.*
+      aeson ==1.5.*
     , attoparsec >=0.13.0.0 && <0.15
-    , base >=4.7 && <4.15
+    , base >=4.7 && <4.16
     , containers >=0.6.0.0 && <0.7
     , formatting >=7.1.0 && <7.2
     , hmatrix >=0.20.0 && <0.21
     , massiv >=0.6.0.0 && <0.7
+    , psqueues >=0.2.7.0 && <0.3
     , rio >=0.1.13.0 && <0.2
   default-language: Haskell2010
 
@@ -104,16 +104,16 @@
   ghc-options: -Wall -Wno-unused-top-binds -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wpartial-fields -Wredundant-constraints -threaded -rtsopts -with-rtsopts=-N
   build-depends:
       ConClusion
-    , PSQueue >=1.1.0.1 && <1.2
     , aeson ==1.5.*
     , attoparsec >=0.13.0.0 && <0.15
-    , base >=4.7 && <4.15
+    , base >=4.7 && <4.16
     , cmdargs >=0.10.0 && <0.11
     , containers >=0.6.0.0 && <0.7
     , formatting >=7.1.0 && <7.2
     , hmatrix >=0.20.0 && <0.21
     , massiv >=0.6.0.0 && <0.7
     , optics >=0.3 && <0.5
+    , psqueues >=0.2.7.0 && <0.3
     , rio >=0.1.13.0 && <0.2
     , text >=1.2.0.0 && <1.3
   default-language: Haskell2010
diff --git a/app/ConClusion.hs b/app/ConClusion.hs
--- a/app/ConClusion.hs
+++ b/app/ConClusion.hs
@@ -257,7 +257,7 @@
           &= help
             "Cluster algorithm to cluster conformeres.\n\
             \  - \"dbscan\" DBScan algorithm. Depends on \n\
-            \             \"measure\", \"distance\" and \"minSize\"\n\
+            \             \"measure\", \"distance\" and \"minsize\"\n\
             \  - \"hca\"    Hierarchical Cluster Analysis. Depends on\n\
             \             \"measure\", \"joinstrat\" and \"distance\"\n\
             \Default: dbscan"
diff --git a/src/ConClusion/Numeric/Statistics.hs b/src/ConClusion/Numeric/Statistics.hs
--- a/src/ConClusion/Numeric/Statistics.hs
+++ b/src/ConClusion/Numeric/Statistics.hs
@@ -41,10 +41,10 @@
 import ConClusion.Numeric.Data hiding (normalise)
 import Data.Aeson hiding (Array)
 import Data.Complex
+import qualified Data.HashPSQ as PQ
 import qualified Data.IntSet as IntSet
 import Data.Massiv.Array as Massiv
 import Data.Massiv.Array.Unsafe as Massiv
-import qualified Data.PSQueue as PQ
 import qualified Numeric.LinearAlgebra as LA
 import RIO hiding (Vector)
 import System.IO.Unsafe (unsafePerformIO)
@@ -103,6 +103,14 @@
     ixVec = makeArrayLinear @D Seq (Sz n') id
     get acc = compute @U . Massiv.map acc
 
+-- | Adjust function for priority queues. Updates the priority at a given key if present.
+pqAdjust :: (Ord k, Hashable k, Ord p) => (p -> p) -> k -> PQ.HashPSQ k p v -> PQ.HashPSQ k p v
+pqAdjust f k q = snd $ PQ.alter f' k q
+  where
+    f' = \op -> case op of
+      Nothing -> (False, Nothing)
+      Just (p, v) -> (False, Just (f p, v))
+
 ----------------------------------------------------------------------------------------------------
 -- Principal Component Analysis
 
@@ -589,7 +597,7 @@
     nNghbr <- nearestNeighbours distMat
 
     let -- Initial priority queue of points. Has the minimum distance of all points.
-        pq = PQ.fromList . Massiv.toList . Massiv.imap (\k (d, _) -> k PQ.:-> d) $ nNghbr
+        pq = PQ.fromList . Massiv.toList . Massiv.imap (\k (d, n) -> (k, d, n)) $ nNghbr
         -- Set of points not joined yet. Initially all points.
         s = IntSet.fromDistinctAscList [0 .. nPoints - 1]
         -- Initial dendrogram accumulator. The vector of all points as their own cluster.
@@ -629,7 +637,7 @@
   -- | List of nearest neighbours for each point.
   MArray (PrimState m) r Ix1 (e, Ix1) ->
   -- | Priority queue with the distances as priorities and the cluster index as keys.
-  PQ.PSQ Ix1 e ->
+  PQ.HashPSQ Ix1 e Ix1 ->
   -- | A set \(S\), that keeps track which clusters have already been joined.
   IntSet ->
   -- | Accumulator of the dendrogram. Should collapse to a singleton vector.
@@ -657,8 +665,11 @@
     -- Redirect neighbours to b, if they previously pointed to a.
     nNghbrU2 <- redirectNeighbours a b newS newDistMat nNghbrU1
 
+    -- Preserve a lower bound in priority queue and update the nearest neighbour list.
+    (nNghbrU3, pqU3) <- updateWithNewBDists b newS newDistMat nNghbrU2 pqU2
+
     -- Update the neighbourlist and priority queue with the new distances to b.
-    (newNNghbr, newPQ) <- updateBNeighbour b s newDistMat nNghbrU2 pqU2
+    (newNNghbr, newPQ) <- updateBNeighbour b s newDistMat nNghbrU3 pqU3
 
     -- If the problem has been reduced to a single cluster the algorithm is done and the final
     -- dendrogram can be obtained from the accumulator at index b. Otherwise join further.
@@ -676,10 +687,10 @@
     Ord e
   ) =>
   MArray (PrimState m) r Ix1 (e, Ix1) ->
-  PQ.PSQ Ix1 e ->
+  PQ.HashPSQ Ix1 e Ix1 ->
   m (Ix1, Ix1, e)
 getJoinCandidates nNghbr pq = do
-  (a PQ.:-> d) <- case PQ.findMin pq of
+  (a, d, _) <- case PQ.findMin pq of
     Nothing -> throwM $ IndexException "Empty priority queue"
     Just v -> return v
   (_, b) <- nNghbr `readM` a
@@ -704,8 +715,8 @@
   IntSet ->
   MArray (PrimState m) r Ix2 e ->
   MArray (PrimState m) r Ix1 (e, Ix1) ->
-  PQ.PSQ Ix1 e ->
-  m (Ix1, Ix1, e, MArray (PrimState m) r Ix1 (e, Ix1), PQ.PSQ Ix1 e)
+  PQ.HashPSQ Ix1 e Ix1 ->
+  m (Ix1, Ix1, e, MArray (PrimState m) r Ix1 (e, Ix1), PQ.HashPSQ Ix1 e Ix1)
 recalculateNghbr (cA, cB, d) s distMat nNghbr pq = do
   dAB <- distMat `readM` (cA :. cB)
   if d == dAB
@@ -718,10 +729,10 @@
       writeM nNghbr cA newNeighbourA
 
       -- Update the priority queue at key cA with the new distance.
-      let newPQ = PQ.adjust (const minDistA) cA pq
+      let newPQ = pqAdjust (const minDistA) cA pq
 
       -- Determine new a, b and d from the updated neighbour list and priority queue.
-      (a PQ.:-> newD) <- case PQ.findMin newPQ of
+      (a, newD, _) <- case PQ.findMin newPQ of
         Nothing -> throwM $ IndexException "Empty priority queue"
         Just v -> return v
       (_, b) <- nNghbr `readM` a
@@ -740,9 +751,9 @@
   Ix1 ->
   e ->
   IntSet ->
-  PQ.PSQ Ix1 e ->
+  PQ.HashPSQ Ix1 e Ix1 ->
   DendroAccM m e ->
-  m (IntSet, PQ.PSQ Ix1 e, DendroAccM m e)
+  m (IntSet, PQ.HashPSQ Ix1 e Ix1, DendroAccM m e)
 joinClusters a b d s pq acc = do
   clA <- acc `readM` a
   let newPQ = PQ.deleteMin pq
@@ -833,7 +844,7 @@
 
 -- | Updates the list of nearest neighbours for all combinations that might have changed by
 -- recalculation with the joined cluster AB at index b.
--- L
+-- L 33-38
 {-# SCC updateWithNewBDists #-}
 updateWithNewBDists ::
   ( MonadThrow m,
@@ -847,8 +858,8 @@
   IntSet ->
   MArray (PrimState m) r Ix2 e ->
   MArray (PrimState m) r Ix1 (e, Ix1) ->
-  PQ.PSQ Ix1 e ->
-  m (MArray (PrimState m) r Ix1 (e, Ix1), PQ.PSQ Ix1 e)
+  PQ.HashPSQ Ix1 e Ix1 ->
+  m (MArray (PrimState m) r Ix1 (e, Ix1), PQ.HashPSQ Ix1 e Ix1)
 updateWithNewBDists b s distMat nNghbr pq = do
   pqT <- newTVarIO pq
   forIO_ ixV $ \ix -> do
@@ -856,11 +867,11 @@
     currentPQ <- readTVarIO pqT
     minDistX <- case PQ.lookup ix currentPQ of
       Nothing -> throwM $ IndexException "Empty priority queue."
-      Just v -> return v
+      Just (p, _v) -> return p
     if dBX < minDistX
       then do
         writeM nNghbr ix (dBX, b)
-        atomically . writeTVar pqT . PQ.adjust (const dBX) ix $ currentPQ
+        atomically . writeTVar pqT . pqAdjust (const dBX) ix $ currentPQ
       else atomically . writeTVar pqT $ currentPQ
 
   newPQ <- readTVarIO pqT
@@ -883,8 +894,8 @@
   IntSet ->
   MArray (PrimState m) r Ix2 e ->
   MArray (PrimState m) r Ix1 (e, Ix1) ->
-  PQ.PSQ Ix1 e ->
-  m (MArray (PrimState m) r Ix1 (e, Ix1), PQ.PSQ Ix1 e)
+  PQ.HashPSQ Ix1 e Ix1 ->
+  m (MArray (PrimState m) r Ix1 (e, Ix1), PQ.HashPSQ Ix1 e Ix1)
 updateBNeighbour b s distMat nNghbr pq =
   if b >= nNeighbours
     then return (nNghbr, pq)
@@ -892,7 +903,7 @@
       rowAB <- searchRow b s distMat >>= unsafeFreeze Par
       newNeighbourB@(distB, neighbourB) <- minimumM rowAB
       writeM nNghbr b newNeighbourB
-      let newPQ = PQ.adjust (const distB) neighbourB pq
+      let newPQ = pqAdjust (const distB) neighbourB pq
       return (nNghbr, newPQ)
   where
     Sz nNeighbours = msize nNghbr
@@ -941,7 +952,8 @@
 searchRow x s dm =
   makeMArray Par (size ixV) $ \ix -> do
     dmIx <- ixV !? ix
-    (dm `readM` (x :. dmIx)) >>= \dist -> return (dist, dmIx)
+    val <- (dm `readM` (x :. dmIx)) >>= \dist -> return (dist, dmIx)
+    return val
   where
     ixV :: Vector U Ix1
     ixV = compute @U . sfilter (> x) . Massiv.fromList @U Par . IntSet.toAscList $ s