diff --git a/ConClusion.cabal b/ConClusion.cabal
--- a/ConClusion.cabal
+++ b/ConClusion.cabal
@@ -1,11 +1,11 @@
 cabal-version: 1.12
 
--- This file has been generated from package.yaml by hpack version 0.34.4.
+-- This file has been generated from package.yaml by hpack version 0.34.6.
 --
 -- see: https://github.com/sol/hpack
 
 name:           ConClusion
-version:        0.1.0
+version:        0.2.0
 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
@@ -27,8 +27,10 @@
 
 library
   exposed-modules:
+      ConClusion.Array.Conversion
+      ConClusion.Array.Util
+      ConClusion.BinaryTree
       ConClusion.Chemistry.Topology
-      ConClusion.Numeric.Data
       ConClusion.Numeric.Statistics
   other-modules:
       Paths_ConClusion
@@ -60,9 +62,9 @@
       NamedFieldPuns
   ghc-options: -Wall -Wno-unused-top-binds -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wpartial-fields -Wredundant-constraints
   build-depends:
-      aeson ==1.5.*
+      aeson >=1.5 && <2.1
     , attoparsec >=0.13.0.0 && <0.15
-    , base >=4.7 && <4.16
+    , base >=4.7 && <4.17
     , containers >=0.6.0.0 && <0.7
     , formatting >=7.1.0 && <7.2
     , hmatrix >=0.20.0 && <0.21
@@ -104,9 +106,9 @@
   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
-    , aeson ==1.5.*
+    , aeson >=1.5 && <2.1
     , attoparsec >=0.13.0.0 && <0.15
-    , base >=4.7 && <4.16
+    , base >=4.7 && <4.17
     , cmdargs >=0.10.0 && <0.11
     , containers >=0.6.0.0 && <0.7
     , formatting >=7.1.0 && <7.2
@@ -115,5 +117,5 @@
     , 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
+    , text >=1.2.0.0 && <2.1
   default-language: Haskell2010
diff --git a/src/ConClusion/Array/Conversion.hs b/src/ConClusion/Array/Conversion.hs
new file mode 100644
--- /dev/null
+++ b/src/ConClusion/Array/Conversion.hs
@@ -0,0 +1,46 @@
+-- |
+-- Module      : ConClusion.Array.Conversion
+-- Description : Type castings and conversions of array types
+-- Copyright   : Phillip Seeber, 2022
+-- License     : AGPL-3
+-- Maintainer  : phillip.seeber@googlemail.com
+-- Stability   : experimental
+-- Portability : POSIX, Windows
+module ConClusion.Array.Conversion
+  ( vecH2M,
+    vecM2H,
+    matH2M,
+    matM2H,
+  )
+where
+
+import Data.Massiv.Array as Massiv hiding (IndexException)
+import Data.Massiv.Array.Manifest.Vector as Massiv
+import Numeric.LinearAlgebra as LA hiding (magnitude, (<>))
+import RIO
+import qualified RIO.Vector.Storable as VectorS
+
+-- | Converts a vector from the HMatrix package to the Massiv representation.
+{-# SCC vecH2M #-}
+vecH2M :: (Element e, Mutable r e, Load r Ix1 e) => VectorS.Vector e -> Massiv.Vector r e
+vecH2M hVec = fromVector' Seq (Sz $ VectorS.length hVec) hVec
+
+-- | Converts a vector from the Massiv representation to the HMatrix representation.
+{-# SCC vecM2H #-}
+vecM2H :: (Manifest r e, Load r Ix1 e, Element e) => Massiv.Vector r e -> LA.Vector e
+vecM2H = Massiv.toVector
+
+-- | Converts a matrix from the HMatrix representation to the Massiv representation.
+{-# SCC matH2M #-}
+matH2M :: (Mutable r e, Load r Ix1 e, Element e) => LA.Matrix e -> Massiv.Matrix r e
+matH2M hMat = Massiv.resize' (Sz $ nRows :. nCols) . vecH2M . LA.flatten $ hMat
+  where
+    nRows = LA.rows hMat
+    nCols = LA.cols hMat
+
+-- | Converts a matrix from Massiv to HMatrix representation.
+{-# SCC matM2H #-}
+matM2H :: (Element e, Manifest r e, Load r Ix1 e) => Massiv.Matrix r e -> LA.Matrix e
+matM2H mMat = LA.reshape nCols . vecM2H . Massiv.flatten $ mMat
+  where
+    Sz (_nRows :. nCols) = Massiv.size mMat
diff --git a/src/ConClusion/Array/Util.hs b/src/ConClusion/Array/Util.hs
new file mode 100644
--- /dev/null
+++ b/src/ConClusion/Array/Util.hs
@@ -0,0 +1,102 @@
+-- |
+-- Module      : ConClusion.Array.Util
+-- Description : Additional tools to work with numerical arrays
+-- Copyright   : Phillip Seeber, 2022
+-- License     : AGPL-3
+-- Maintainer  : phillip.seeber@googlemail.com
+-- Stability   : experimental
+-- Portability : POSIX, Windows
+module ConClusion.Array.Util
+  ( IndexException (..),
+    magnitude,
+    normalise,
+    angle,
+    minDistAt,
+    minDistAtVec,
+    iMinimumM,
+  )
+where
+
+import Data.Massiv.Array as Massiv hiding (IndexException)
+import RIO
+import System.IO.Unsafe (unsafePerformIO)
+
+-- | Exception regarding indexing in some kind of aaray.
+newtype IndexException = IndexException String deriving (Show)
+
+instance Exception IndexException
+
+-- | Magnitude of a vector (length).
+magnitude :: (Massiv.Numeric r e, Source r e, Floating e) => Massiv.Vector r e -> e
+magnitude v = sqrt $ v !.! v
+
+-- | Normalise a vector.
+normalise :: (Massiv.Numeric r e, Source r e, Floating e) => Massiv.Vector r e -> Massiv.Vector r e
+normalise v = v .* (1 / magnitude v)
+
+-- | Angle between two vectors.
+angle :: (Massiv.Numeric r e, Source r e, Floating e) => Massiv.Vector r e -> Massiv.Vector r e -> e
+angle a b = acos $ a !.! b / (magnitude a * magnitude b)
+
+-- | Find the minimal distance in a distance matrix, which is not the main diagonal.
+{-# SCC minDistAt #-}
+minDistAt ::
+  ( Manifest r e,
+    MonadThrow m,
+    Ord e
+  ) =>
+  Massiv.Matrix r e ->
+  m (e, Ix2)
+minDistAt arr
+  | isEmpty arr = throwM $ SizeEmptyException (Massiv.size arr)
+  | otherwise = return . unsafePerformIO $ ifoldlP minFold start chFold start arr
+  where
+    ix0 = pureIndex 0
+    e0 = arr Massiv.! ix0
+    start = (e0, ix0)
+    minFold acc@(eA, _) ix@(m :. n) e = if e < eA && m > n then (e, ix) else acc
+    chFold acc@(eA, _) ch@(e, _) = if e < eA then ch else acc
+
+-- | Find the minimal element of a vector, which is at a larger than the supplied index.
+minDistAtVec ::
+  ( Manifest r e,
+    MonadThrow m,
+    Ord e
+  ) =>
+  Ix1 ->
+  Massiv.Vector r e ->
+  m (e, Ix1)
+minDistAtVec ixStart vec
+  | isEmpty vec = throwM $ SizeEmptyException (Massiv.size vec)
+  | ixStart >= nElems = throwM $ IndexOutOfBoundsException (Sz nElems) ixStart
+  | otherwise = do
+    let (minE, minIx) = unsafePerformIO $ ifoldlP minFold startAcc chFold startAcc searchVec
+    return (minE, minIx + ixStart + 1)
+  where
+    Sz nElems = Massiv.size vec
+    searchVec = Massiv.drop (Sz $ ixStart + 1) vec
+    ix0 = 0
+    e0 = searchVec Massiv.! ix0
+    startAcc = (e0, ix0)
+    minFold acc@(eA, _) ix e = if e < eA then (e, ix) else acc
+    chFold acc ch = min acc ch
+
+-- | Like 'Massiv.minimumM' but also returns the index of the minimal element.
+iMinimumM ::
+  ( Manifest r a,
+    MonadThrow m,
+    Index ix,
+    Ord a
+  ) =>
+  Array r ix a ->
+  m (a, ix)
+iMinimumM arr
+  | isEmpty arr = throwM $ SizeEmptyException (Massiv.size arr)
+  | otherwise = return . unsafePerformIO $ ifoldlP minFold start chFold start arr
+  where
+    ix0 = pureIndex 0
+    e0 = arr Massiv.! ix0
+    start = (e0, ix0)
+
+    minFold acc@(eA, _) ix e = if e < eA then (e, ix) else acc
+    chFold acc@(eA, _) ch@(e, _) = if e < eA then ch else acc
diff --git a/src/ConClusion/BinaryTree.hs b/src/ConClusion/BinaryTree.hs
new file mode 100644
--- /dev/null
+++ b/src/ConClusion/BinaryTree.hs
@@ -0,0 +1,60 @@
+-- |
+-- Module      : ConClusion.BinaryTree
+-- Description : Custom binary tree type with some special functions
+-- Copyright   : Phillip Seeber, 2022
+-- License     : AGPL-3
+-- Maintainer  : phillip.seeber@googlemail.com
+-- Stability   : experimental
+-- Portability : POSIX, Windows
+module ConClusion.BinaryTree
+  ( BinTree (..),
+    root,
+    takeBranchesWhile,
+    takeLeafyBranchesWhile,
+  )
+where
+
+import Data.Aeson hiding (Array)
+import Data.Massiv.Array as Massiv hiding (IndexException)
+import RIO
+
+-- | A binary tree.
+data BinTree e = Leaf e | Node e (BinTree e) (BinTree e)
+  deriving (Eq, Show, Generic)
+
+instance (FromJSON e) => FromJSON (BinTree e)
+
+instance (ToJSON e) => ToJSON (BinTree e)
+
+instance Functor BinTree where
+  fmap f (Leaf a) = Leaf (f a)
+  fmap f (Node a l r) = Node (f a) (fmap f l) (fmap f r)
+
+-- | Look at the root of a binary tree.
+root :: BinTree e -> e
+root (Leaf e) = e
+root (Node e _ _) = e
+
+-- | Steps down each branch of a tree until some criterion is satisfied or the end of the branch is
+-- reached. Each end of the branch is added to a result.
+takeBranchesWhile :: (a -> Bool) -> BinTree a -> Massiv.Vector DL a
+takeBranchesWhile chk tree = go tree (Massiv.empty @DL)
+  where
+    go (Leaf v) acc = if chk v then acc `snoc` v else acc
+    go (Node v l r) acc =
+      let vAcc = if chk v then acc `snoc` v else acc
+          lAcc = go l vAcc
+          rAcc = go r lAcc
+       in if chk v then rAcc else vAcc
+
+-- | Takes the first value in each branch, that does not fullfill the criterion anymore and adds it
+-- to the result. Terminal leafes of the branches are always taken.
+takeLeafyBranchesWhile :: (a -> Bool) -> BinTree a -> Massiv.Vector DL a
+takeLeafyBranchesWhile chk tree = go tree (Massiv.empty @DL)
+  where
+    go (Leaf v) acc = acc `snoc` v
+    go (Node v l r) acc =
+      let vAcc = if chk v then acc else acc `snoc` v
+          lAcc = go l vAcc
+          rAcc = go r lAcc
+       in if chk v then rAcc else vAcc
diff --git a/src/ConClusion/Chemistry/Topology.hs b/src/ConClusion/Chemistry/Topology.hs
--- a/src/ConClusion/Chemistry/Topology.hs
+++ b/src/ConClusion/Chemistry/Topology.hs
@@ -26,8 +26,9 @@
   )
 where
 
-import ConClusion.Numeric.Data hiding (angle)
-import qualified ConClusion.Numeric.Data as ND
+import ConClusion.Array.Conversion
+import ConClusion.Array.Util hiding (angle)
+import qualified ConClusion.Array.Util as ArrayUtil
 import Data.Attoparsec.Text hiding (D)
 import Data.Foldable
 import Data.Massiv.Array as Massiv hiding (B, D, forM)
@@ -135,7 +136,7 @@
     coordC <- compute @U <$> (coordinates !?> c)
     vecAB <- coordA .-. coordB
     vecCB <- coordC .-. coordB
-    return $ ND.angle vecAB vecCB
+    return $ ArrayUtil.angle vecAB vecCB
 
 -- | Calculates the dihedral angle defined by four atoms. Respects rotation direction. Obtains the
 -- result in radian.
@@ -157,7 +158,7 @@
           if vecH2M normVecRot !.! vecAB < 0
             then -1
             else 1
-    return $ rotDir * ND.angle planeABC planeBCD
+    return $ rotDir * ArrayUtil.angle planeABC planeBCD
 
 -- | Calculates a metric value of the dihedral angle defined by four atoms. This must create 2
 -- values in the feature matrix, instead of one.
diff --git a/src/ConClusion/Numeric/Data.hs b/src/ConClusion/Numeric/Data.hs
deleted file mode 100644
--- a/src/ConClusion/Numeric/Data.hs
+++ /dev/null
@@ -1,198 +0,0 @@
--- |
--- Module      : ConClusion.Numeric.Data
--- Description : Type castings, conversions and utilities
--- Copyright   : Phillip Seeber, 2021
--- License     : AGPL-3
--- Maintainer  : phillip.seeber@googlemail.com
--- Stability   : experimental
--- Portability : POSIX, Windows
-module ConClusion.Numeric.Data
-  ( -- * Conversion of array types.
-
-    -- Conversion between HMatrix and Massiv's array types
-    IndexException (..),
-    vecH2M,
-    vecM2H,
-    matH2M,
-    matM2H,
-
-    -- * Array Processing
-    magnitude,
-    normalise,
-    angle,
-    minDistAt,
-    minDistAtVec,
-    iMinimumM,
-
-    -- * Utilities
-    printMat,
-
-    -- * Binary Trees
-    BinTree (..),
-    root,
-    takeBranchesWhile,
-    takeLeafyBranchesWhile,
-  )
-where
-
-import Data.Aeson hiding (Array)
-import Data.Massiv.Array as Massiv hiding (IndexException)
-import Data.Massiv.Array.Manifest.Vector as Massiv
-import Formatting
-import Numeric.LinearAlgebra as LA hiding (magnitude, (<>))
-import RIO
-import qualified RIO.Vector.Storable as VectorS
-import System.IO.Unsafe (unsafePerformIO)
-
--- | Exception regarding indexing in some kind of aaray.
-newtype IndexException = IndexException String deriving (Show)
-
-instance Exception IndexException
-
--- | Converts a vector from the HMatrix package to the Massiv representation.
-{-# SCC vecH2M #-}
-vecH2M :: (Element e, Mutable r e, Load r Ix1 e) => VectorS.Vector e -> Massiv.Vector r e
-vecH2M hVec = fromVector' Seq (Sz $ VectorS.length hVec) hVec
-
--- | Converts a vector from the Massiv representation to the HMatrix representation.
-{-# SCC vecM2H #-}
-vecM2H :: (Manifest r e, Load r Ix1 e, Element e) => Massiv.Vector r e -> LA.Vector e
-vecM2H = Massiv.toVector
-
--- | Converts a matrix from the HMatrix representation to the Massiv representation.
-{-# SCC matH2M #-}
-matH2M :: (Mutable r e, Load r Ix1 e, Element e) => LA.Matrix e -> Massiv.Matrix r e
-matH2M hMat = Massiv.resize' (Sz $ nRows :. nCols) . vecH2M . LA.flatten $ hMat
-  where
-    nRows = LA.rows hMat
-    nCols = LA.cols hMat
-
--- | Converts a matrix from Massiv to HMatrix representation.
-{-# SCC matM2H #-}
-matM2H :: (Element e, Manifest r e, Load r Ix1 e) => Massiv.Matrix r e -> LA.Matrix e
-matM2H mMat = LA.reshape nCols . vecM2H . Massiv.flatten $ mMat
-  where
-    Sz (_nRows :. nCols) = Massiv.size mMat
-
--- | Magnitude of a vector (length).
-magnitude :: (Massiv.Numeric r e, Source r e, Floating e) => Massiv.Vector r e -> e
-magnitude v = sqrt $ v !.! v
-
--- | Normalise a vector.
-normalise :: (Massiv.Numeric r e, Source r e, Floating e) => Massiv.Vector r e -> Massiv.Vector r e
-normalise v = v .* (1 / magnitude v)
-
--- | Angle between two vectors.
-angle :: (Massiv.Numeric r e, Source r e, Floating e) => Massiv.Vector r e -> Massiv.Vector r e -> e
-angle a b = acos $ a !.! b / (magnitude a * magnitude b)
-
--- | Find the minimal distance in a distance matrix, which is not the main diagonal.
-{-# SCC minDistAt #-}
-minDistAt ::
-  ( Manifest r e,
-    MonadThrow m,
-    Ord e
-  ) =>
-  Massiv.Matrix r e ->
-  m (e, Ix2)
-minDistAt arr
-  | isEmpty arr = throwM $ SizeEmptyException (Massiv.size arr)
-  | otherwise = return . unsafePerformIO $ ifoldlP minFold start chFold start arr
-  where
-    ix0 = pureIndex 0
-    e0 = arr Massiv.! ix0
-    start = (e0, ix0)
-    minFold acc@(eA, _) ix@(m :. n) e = if e < eA && m > n then (e, ix) else acc
-    chFold acc@(eA, _) ch@(e, _) = if e < eA then ch else acc
-
--- | Find the minimal element of a vector, which is at a larger than the supplied index.
-minDistAtVec ::
-  ( Manifest r e,
-    MonadThrow m,
-    Ord e
-  ) =>
-  Ix1 ->
-  Massiv.Vector r e ->
-  m (e, Ix1)
-minDistAtVec ixStart vec
-  | isEmpty vec = throwM $ SizeEmptyException (Massiv.size vec)
-  | ixStart >= nElems = throwM $ IndexOutOfBoundsException (Sz nElems) ixStart
-  | otherwise = do
-    let (minE, minIx) = unsafePerformIO $ ifoldlP minFold startAcc chFold startAcc searchVec
-    return (minE, minIx + ixStart + 1)
-  where
-    Sz nElems = Massiv.size vec
-    searchVec = Massiv.drop (Sz $ ixStart + 1) vec
-    ix0 = 0
-    e0 = searchVec Massiv.! ix0
-    startAcc = (e0, ix0)
-    minFold acc@(eA, _) ix e = if e < eA then (e, ix) else acc
-    chFold acc ch = min acc ch
-
--- | Like 'Massiv.minimumM' but also returns the index of the minimal element.
-iMinimumM ::
-  ( Manifest r a,
-    MonadThrow m,
-    Index ix,
-    Ord a
-  ) =>
-  Array r ix a ->
-  m (a, ix)
-iMinimumM arr
-  | isEmpty arr = throwM $ SizeEmptyException (Massiv.size arr)
-  | otherwise = return . unsafePerformIO $ ifoldlP minFold start chFold start arr
-  where
-    ix0 = pureIndex 0
-    e0 = arr Massiv.! ix0
-    start = (e0, ix0)
-
-    minFold acc@(eA, _) ix e = if e < eA then (e, ix) else acc
-    chFold acc@(eA, _) ch@(e, _) = if e < eA then ch else acc
-
--- | Quickly print a matrix with numerical values
-printMat :: (Source r e, Real e) => Massiv.Matrix r e -> Massiv.Matrix D Text
-printMat mat = Massiv.map (sformat (left 4 ' ' %. fixed 2)) mat
-
-----------------------------------------------------------------------------------------------------
--- Binary Trees.
-
--- | A binary tree.
-data BinTree e = Leaf e | Node e (BinTree e) (BinTree e)
-  deriving (Eq, Show, Generic)
-
-instance (FromJSON e) => FromJSON (BinTree e)
-
-instance (ToJSON e) => ToJSON (BinTree e)
-
-instance Functor BinTree where
-  fmap f (Leaf a) = Leaf (f a)
-  fmap f (Node a l r) = Node (f a) (fmap f l) (fmap f r)
-
--- | Look at the root of a binary tree.
-root :: BinTree e -> e
-root (Leaf e) = e
-root (Node e _ _) = e
-
--- | Steps down each branch of a tree until some criterion is satisfied or the end of the branch is
--- reached. Each end of the branch is added to a result.
-takeBranchesWhile :: (a -> Bool) -> BinTree a -> Massiv.Vector DL a
-takeBranchesWhile chk tree = go tree (Massiv.empty @DL)
-  where
-    go (Leaf v) acc = if chk v then acc `snoc` v else acc
-    go (Node v l r) acc =
-      let vAcc = if chk v then acc `snoc` v else acc
-          lAcc = go l vAcc
-          rAcc = go r lAcc
-       in if chk v then rAcc else vAcc
-
--- | Takes the first value in each branch, that does not fullfill the criterion anymore and adds it
--- to the result. Terminal leafes of the branches are always taken.
-takeLeafyBranchesWhile :: (a -> Bool) -> BinTree a -> Massiv.Vector DL a
-takeLeafyBranchesWhile chk tree = go tree (Massiv.empty @DL)
-  where
-    go (Leaf v) acc = acc `snoc` v
-    go (Node v l r) acc =
-      let vAcc = if chk v then acc else acc `snoc` v
-          lAcc = go l vAcc
-          rAcc = go r lAcc
-       in if chk v then rAcc else vAcc
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
@@ -38,7 +38,9 @@
   )
 where
 
-import ConClusion.Numeric.Data hiding (normalise)
+import ConClusion.Array.Conversion
+import ConClusion.Array.Util hiding (normalise)
+import ConClusion.BinaryTree
 import Data.Aeson hiding (Array)
 import Data.Complex
 import qualified Data.HashPSQ as PQ
@@ -64,8 +66,6 @@
     Load r1 Ix1 (Complex Double),
     Load r2 Ix1 (Complex Double),
     Load r3 Ix1 e,
-    -- Resize r3 Ix2,
-    Load r3 Ix2 e,
     MonadThrow m
   ) =>
   Matrix r3 e ->
@@ -579,9 +579,6 @@
   ( MonadThrow m,
     Mutable r e,
     Mutable r (e, Ix1),
-    -- Manifest (R r) Ix1 e,
-    -- OuterSlice r Ix2 e,
-    Shape r Ix1,
     Load r Ix1 e,
     Ord e,
     Unbox e,