diff --git a/Changelog.md b/Changelog.md
--- a/Changelog.md
+++ b/Changelog.md
@@ -1,5 +1,9 @@
 # Changelog
 
+## 0.1.0
+  - updates to Massiv 1.0.0.0; changes lots of array types.
+  - infrastructure updates
+
 ## 0.0.2
   - switching from `PSQueue` to `psqueues` as the former is not maintained
   - raise upper bounds of base, allows GHC 9.0
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.2
+version:        0.1.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
@@ -66,7 +66,7 @@
     , 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
+    , massiv >=1.0.0.0 && <1.1
     , psqueues >=0.2.7.0 && <0.3
     , rio >=0.1.13.0 && <0.2
   default-language: Haskell2010
@@ -111,7 +111,7 @@
     , 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
+    , massiv >=1.0.0.0 && <1.1
     , optics >=0.3 && <0.5
     , psqueues >=0.2.7.0 && <0.3
     , rio >=0.1.13.0 && <0.2
diff --git a/app/ConClusion.hs b/app/ConClusion.hs
--- a/app/ConClusion.hs
+++ b/app/ConClusion.hs
@@ -621,7 +621,7 @@
   valueLines <- case clusters of
     Nothing -> return . Massiv.fold . Massiv.map (writeLine . Left) . innerSlices $ mat
     Just cl -> do
-      let points = innerSlices mat :: Vector Massiv.D (Vector M Double)
+      let points = innerSlices mat
           labeledCl = compute . Massiv.imap (,) $ cl :: Vector Massiv.B (Int, IntSet)
 
       -- Anotate all points with their assignment to a cluster number. If a point is from DBScan, it
@@ -666,12 +666,12 @@
           <> bformat (left cw ' ' F.%. (builder F.% iForm F.% "," F.% iForm F.% "," F.% iForm F.% "," F.% iForm)) "A cos" a b c d
 
     -- Print a vector of numeric values.
-    doubleVecB :: Source r Ix1 Double => Vector r Double -> TB.Builder
+    doubleVecB :: Source r Double => Vector r Double -> TB.Builder
     doubleVecB v = Massiv.fold . Massiv.map dForm $ v
 
     -- Writer for a value line.
     writeLine ::
-      (Source r Ix1 Double) => Either (Vector r Double) (Maybe Int, Vector r Double) -> TB.Builder
+      (Source r Double) => Either (Vector r Double) (Maybe Int, Vector r Double) -> TB.Builder
     writeLine (Left vec) = doubleVecB vec <> "\n"
     writeLine (Right (cl, vec)) = bformat iForm (fromMaybe (-1) cl) <> doubleVecB vec <> "\n"
 
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
@@ -159,7 +159,7 @@
             else 1
     return $ rotDir * ND.angle planeABC planeBCD
 
--- | Calculates the a metric value of the dihedral angle defined by four atoms. This must create 2
+-- | Calculates a metric value of the dihedral angle defined by four atoms. This must create 2
 -- values in the feature matrix, instead of one.
 -- See <https://onlinelibrary.wiley.com/doi/full/10.1002/prot.20310)>
 dihedral :: MonadThrow m => D -> Molecule -> m (Double, Double)
@@ -207,4 +207,4 @@
   m (Massiv.Matrix DL Double)
 getFeatures sel trj = traverse toCols trj >>= concatM (Dim 1)
   where
-    toCols v = expandInner @Ix2 (Sz 1) const . compute @U <$> getFeature sel v
+    toCols v = expandInner @U @Ix2 (Sz 1) const . compute @U <$> getFeature sel v
diff --git a/src/ConClusion/Numeric/Data.hs b/src/ConClusion/Numeric/Data.hs
--- a/src/ConClusion/Numeric/Data.hs
+++ b/src/ConClusion/Numeric/Data.hs
@@ -51,17 +51,17 @@
 
 -- | Converts a vector from the HMatrix package to the Massiv representation.
 {-# SCC vecH2M #-}
-vecH2M :: (Element e, Mutable r Ix1 e) => VectorS.Vector e -> Massiv.Vector r e
+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 Ix1 e, Element e) => Massiv.Vector r e -> LA.Vector e
+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 Ix1 e, Element e) => LA.Matrix e -> Massiv.Matrix r e
+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
@@ -69,27 +69,27 @@
 
 -- | Converts a matrix from Massiv to HMatrix representation.
 {-# SCC matM2H #-}
-matM2H :: (Manifest r Ix1 e, Element e, Resize r Ix2, Load r Ix2 e) => Massiv.Matrix r e -> LA.Matrix e
+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 Ix1 e, Floating e) => Massiv.Vector r e -> e
+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 Ix1 e, Floating e) => Massiv.Vector r e -> Massiv.Vector r e
+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 Ix1 e, Floating e) => Massiv.Vector r e -> Massiv.Vector r e -> e
+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 Ix2 e,
+  ( Manifest r e,
     MonadThrow m,
     Ord e
   ) =>
@@ -107,7 +107,7 @@
 
 -- | Find the minimal element of a vector, which is at a larger than the supplied index.
 minDistAtVec ::
-  ( Manifest r Ix1 e,
+  ( Manifest r e,
     MonadThrow m,
     Ord e
   ) =>
@@ -131,8 +131,9 @@
 
 -- | Like 'Massiv.minimumM' but also returns the index of the minimal element.
 iMinimumM ::
-  ( Manifest r ix a,
+  ( Manifest r a,
     MonadThrow m,
+    Index ix,
     Ord a
   ) =>
   Array r ix a ->
@@ -149,7 +150,7 @@
     chFold acc@(eA, _) ch@(e, _) = if e < eA then ch else acc
 
 -- | Quickly print a matrix with numerical values
-printMat :: (Source r Ix2 e, Real e) => Massiv.Matrix r e -> Massiv.Matrix D Text
+printMat :: (Source r e, Real e) => Massiv.Matrix r e -> Massiv.Matrix D Text
 printMat mat = Massiv.map (sformat (left 4 ' ' %. fixed 2)) mat
 
 ----------------------------------------------------------------------------------------------------
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
@@ -55,11 +55,16 @@
 -- | Solves eigenvalue problem of a square matrix and obtains its eigenvalues and eigenvectors.
 {-# SCC eig #-}
 eig ::
-  ( Mutable r1 Ix1 (Complex Double),
-    Mutable r2 Ix1 (Complex Double),
+  ( -- Mutable r1 Ix1 (Complex Double),
+    -- Mutable r2 Ix1 (Complex Double),
     LA.Field e,
-    Manifest r3 Ix1 e,
-    Resize r3 Ix2,
+    Manifest r3 e,
+    Manifest r1 (Complex Double),
+    Manifest r2 (Complex Double),
+    Load r1 Ix1 (Complex Double),
+    Load r2 Ix1 (Complex Double),
+    Load r3 Ix1 e,
+    -- Resize r3 Ix2,
     Load r3 Ix2 e,
     MonadThrow m
   ) =>
@@ -78,10 +83,10 @@
 eigSort ::
   ( Load r2 Ix2 e,
     MonadThrow m,
-    Source r1 Ix1 e,
-    Source r2 Ix2 e,
-    Mutable r1 Ix1 e,
-    Mutable r2 Ix2 e,
+    Source r1 e,
+    Source r2 e,
+    Mutable r1 e,
+    Mutable r2 e,
     Unbox e,
     Ord e
   ) =>
@@ -141,9 +146,9 @@
 -- constructed from the eigenvectors associated to the largest eigenvalues.
 {-# SCC transformToPCABasis #-}
 transformToPCABasis ::
-  ( Source (R r) Ix2 e,
-    Extract r Ix2 e,
-    Mutable r Ix2 e,
+  ( -- Source (R r) Ix2 e,
+    -- Extract r Ix2 e,
+    Mutable r e,
     Numeric r e,
     MonadThrow m
   ) =>
@@ -177,10 +182,9 @@
 {-# SCC pca #-}
 pca ::
   ( Numeric r Double,
-    Mutable r Ix2 Double,
-    Manifest r Ix1 Double,
-    Source (R r) Ix2 Double,
-    Extract r Ix2 Double,
+    Mutable r Double,
+    Load r Ix1 Double,
+    Load r Ix2 Double,
     MonadThrow m
   ) =>
   -- | Dimensionalty after PCA transformation.
@@ -238,11 +242,11 @@
 -- mean deviation form.
 {-# SCC meanDeviation #-}
 meanDeviation ::
-  ( Source r Ix2 e,
+  ( Source r e,
     Fractional e,
     Unbox e,
     Numeric r e,
-    Mutable r Ix2 e
+    Mutable r e
   ) =>
   Matrix r e ->
   Matrix r e
@@ -260,7 +264,7 @@
 --
 -- The feature matrix should be in mean deviation form, see 'meanDeviation'.
 {-# SCC covariance #-}
-covariance :: (Numeric r e, Mutable r Ix2 e, Fractional e) => Matrix r e -> Matrix r e
+covariance :: (Numeric r e, Mutable r e, Fractional e) => Matrix r e -> Matrix r e
 covariance x = (1 / (fromIntegral n - 1)) *. (x !><! (compute . transpose $ x))
   where
     Sz (_ :. n) = size x
@@ -271,15 +275,14 @@
     Unbox e,
     Numeric r e,
     Fractional e,
-    Source r Ix2 e,
-    Mutable r Ix2 e
+    Mutable r e
   ) =>
   Array r Ix2 e ->
   Array r Ix2 e
 normalise mat =
   let absMat = Massiv.map abs mat
       maxPerRow = compute @U . foldlInner max 0 $ absMat
-      divMat = compute . Massiv.map (1 /) . expandInner @Ix2 (Sz n) const $ maxPerRow
+      divMat = compute . Massiv.map (1 /) . expandInner @U @Ix2 (Sz n) const $ maxPerRow
    in divMat !*! mat
   where
     Sz (_ :. n) = size mat
@@ -292,7 +295,7 @@
 
 -- | Builds the distance measures in a permutation matrix/distance matrix.
 buildDistMat ::
-  (Mutable r Ix2 e) =>
+  (Manifest r e) =>
   -- | Zip function to combine the elements of vectors \(\mathbf{a}\) \(\mathbf{b}\). Usually @(-)@.
   -- \( f(\mathbf{a}_i, \mathbf{b}_i) = \mathbf{c} \)
   (e -> e -> a) ->
@@ -305,7 +308,7 @@
   -- | Resulting distance matrix.
   Matrix D a
 buildDistMat zipFn foldFn acc mat =
-  let a = transposeOuter . expandOuter @Ix3 (Sz n) const $ mat
+  let a = transposeOuter @D @Ix3 . expandOuter (Sz n) const $ mat
       b = transposeInner a
       ab = Massiv.zipWith zipFn a b
       d = foldlInner foldFn acc ab
@@ -318,7 +321,7 @@
 --   d(\mathbf{a}, \mathbf{b}) = \left( \sum \limits_{i=1}^n \lvert \mathbf{a}_i - \mathbf{b}_i \rvert ^p \right) ^ \frac{1}{p}
 -- \]
 {-# SCC lpNorm #-}
-lpNorm :: (Mutable r Ix2 e, Floating e) => Int -> DistFn r e
+lpNorm :: (Mutable r e, Floating e) => Int -> DistFn r e
 lpNorm p = compute . buildDistMat zipFn foldFn 0
   where
     zipFn a b = (^ p) . abs $ a - b
@@ -329,7 +332,7 @@
 --   d(\mathbf{a}, \mathbf{b}) = \sum \limits_{i=1}^n \lvert \mathbf{a}_i - \mathbf{b}_i \rvert
 -- \]
 {-# SCC manhattan #-}
-manhattan :: (Mutable r Ix2 e, Floating e) => DistFn r e
+manhattan :: (Mutable r e, Floating e) => DistFn r e
 manhattan = lpNorm 1
 
 -- | The Euclidean distance between two vectors. Specialisation of the \(L_p\) norm for \(p = 2\).
@@ -337,7 +340,7 @@
 --   d(\mathbf{a}, \mathbf{b}) = \sqrt{\sum \limits_{i=1}^n (\mathbf{a}_i - \mathbf{b}_i)^2}
 -- \]
 {-# SCC euclidean #-}
-euclidean :: (Mutable r Ix2 e, Floating e) => DistFn r e
+euclidean :: (Mutable r e, Floating e) => DistFn r e
 euclidean = lpNorm 2
 
 -- | Mahalanobis distance between points. Suitable for non correlated axes.
@@ -346,9 +349,9 @@
 -- \]
 -- where \(\mathbf{S}\) is the covariance matrix.
 {-# SCC mahalanobis #-}
-mahalanobis :: (Unbox e, Numeric r e, Mutable r Ix2 e, Mutable r Ix1 e, LA.Field e) => DistFn r e
+mahalanobis :: (Unbox e, Numeric r e, Mutable r e, LA.Field e, Load r Ix1 e) => DistFn r e
 mahalanobis points =
-  let a = transposeOuter . expandOuter @Ix3 (Sz n) const $ points
+  let a = transposeOuter @D @Ix3 . expandOuter (Sz n) const $ points
       b = transposeInner a
       abDiff = compute @U $ a !-! b
       ixArray = makeArray @U @Ix2 @Ix2 Par (Sz $ n :. n) id
@@ -384,7 +387,7 @@
     Num e,
     Typeable e,
     Show e,
-    Source r Ix2 e
+    Source r e
   ) =>
   -- | Distance measure to build the distance matrix of all points.
   DistFn r e ->
@@ -416,7 +419,7 @@
     -- Construct the overlap matrix of all clusters.
     compareSets :: (IntSet -> IntSet -> Bool) -> Vector B IntSet -> Matrix D Bool
     compareSets fn clVec =
-      let a = expandOuter @Ix2 sz const clVec
+      let a = expandOuter sz const clVec
           b = transpose a
           compareMat = Massiv.zipWith fn a b
        in compareMat
@@ -574,11 +577,12 @@
 {-# SCC hca #-}
 hca ::
   ( MonadThrow m,
-    Mutable r Ix1 e,
-    Mutable r Ix2 e,
-    Mutable r Ix1 (e, Ix1),
-    Manifest (R r) Ix1 e,
-    OuterSlice r Ix2 e,
+    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,
     Fractional e
@@ -588,7 +592,7 @@
   Matrix r e ->
   m (Dendrogram e)
 hca distFn joinStrat points
-  | Massiv.isEmpty points = throw $ SizeEmptyException (Sz nPoints)
+  | Massiv.isEmpty points = throwM $ SizeEmptyException (Sz nPoints)
   | otherwise = do
     let -- The distance matrix from the points.
         distMat = distFn points
@@ -623,10 +627,11 @@
     PrimMonad m,
     MonadUnliftIO m,
     PrimState m ~ RealWorld,
-    Mutable r Ix2 e,
-    OuterSlice r Ix2 e,
-    Manifest (R r) Ix1 e,
-    Mutable r Ix1 (e, Ix1),
+    Mutable r e,
+    -- OuterSlice r Ix2 e,
+    -- Manifest (R r) Ix1 e,
+    Mutable r (e, Ix1),
+    Shape r Ix1,
     Fractional e,
     Ord e
   ) =>
@@ -645,7 +650,7 @@
   -- | The final dendrogram, after all clusters have been joined.
   m (Dendrogram e)
 agglomerate joinStrat distMat nNghbr pq s dendroAcc
-  | IntSet.null s = throw $ IndexException "No clusters left. This must never happen."
+  | IntSet.null s = throwM $ IndexException "No clusters left. This must never happen."
   | otherwise = do
     -- Obtain candidates for the two clusters to join and the minimal distance in the priority queue.
     candidates <- getJoinCandidates nNghbr pq
@@ -683,7 +688,7 @@
 getJoinCandidates ::
   ( MonadThrow m,
     PrimMonad m,
-    Mutable r Ix1 (e, Ix1),
+    Mutable r (e, Ix1),
     Ord e
   ) =>
   MArray (PrimState m) r Ix1 (e, Ix1) ->
@@ -705,10 +710,11 @@
     PrimMonad m,
     MonadUnliftIO m,
     PrimState m ~ RealWorld,
-    OuterSlice r Ix2 e,
-    Manifest (R r) Ix1 e,
-    Mutable r Ix1 (e, Ix1),
-    Mutable r Ix2 e,
+    -- OuterSlice r Ix2 e,
+    -- Manifest (R r) Ix1 e,
+    Mutable r (e, Ix1),
+    Mutable r e,
+    Shape r Ix1,
     Ord e
   ) =>
   (Ix1, Ix1, e) ->
@@ -782,7 +788,7 @@
   ( MonadThrow m,
     PrimMonad m,
     MonadUnliftIO m,
-    Mutable r Ix2 e,
+    Mutable r e,
     Fractional e
   ) =>
   JoinStrat e ->
@@ -806,8 +812,8 @@
       modifyM_ distMat (\dBX -> return $ lanceWilliams js nA nB nX dAB dAX dBX) (b :. ix)
     return distMat
   where
-    Sz (mDM :. nDM) = msize distMat
-    Sz nCl = msize dendroAcc
+    Sz (mDM :. nDM) = sizeOfMArray distMat
+    Sz nCl = sizeOfMArray dendroAcc
     ixV = Massiv.fromList @U Par . IntSet.toAscList . IntSet.delete b $ s
     clSize i = IntSet.size . cluster . root . unDendro <$> dendroAcc `readM` i
 
@@ -819,8 +825,8 @@
   ( MonadThrow m,
     PrimMonad m,
     MonadUnliftIO m,
-    Mutable r Ix1 (e, Ix1),
-    Mutable r Ix2 e
+    Mutable r (e, Ix1),
+    Mutable r e
   ) =>
   Ix1 ->
   Ix1 ->
@@ -850,8 +856,8 @@
   ( MonadThrow m,
     MonadUnliftIO m,
     PrimMonad m,
-    Mutable r Ix2 e,
-    Mutable r Ix1 (e, Ix1),
+    Mutable r e,
+    Mutable r (e, Ix1),
     Ord e
   ) =>
   Ix1 ->
@@ -879,15 +885,17 @@
   where
     ixV = compute @U . Massiv.sfilter (< b) . Massiv.fromList @U Par . IntSet.toAscList $ s
 
---  | Updates the list of nearest neighbours and the priority queue at key b.
+-- | Updates the list of nearest neighbours and the priority queue at key b.
 -- L 39-40
 {-# SCC updateBNeighbour #-}
 updateBNeighbour ::
   ( MonadThrow m,
     PrimMonad m,
+    RealWorld ~ PrimState m,
     MonadUnliftIO m,
-    Mutable r Ix1 (e, Ix1),
-    Mutable r Ix2 e,
+    Mutable r (e, Ix1),
+    Mutable r e,
+    Shape r Ix1,
     Ord e
   ) =>
   Ix1 ->
@@ -906,7 +914,7 @@
       let newPQ = pqAdjust (const distB) neighbourB pq
       return (nNghbr, newPQ)
   where
-    Sz nNeighbours = msize nNghbr
+    Sz nNeighbours = sizeOfMArray nNghbr
 
 -- | Find the nearest neighbour for each point from a distance matrix. For each point it stores the
 -- minimum distance and the index of the other point, that is the nearest neighbour but at a higher
@@ -914,10 +922,11 @@
 {-# SCC nearestNeighbours #-}
 nearestNeighbours ::
   ( MonadThrow m,
-    Mutable r Ix1 e,
-    Mutable r Ix1 (e, Ix1),
-    OuterSlice r Ix2 e,
-    Source (R r) Ix1 e,
+    Mutable r e,
+    Mutable r (e, Ix1),
+    Load r Ix1 e,
+    -- OuterSlice r Ix2 e,
+    -- Source (R r) Ix1 e,
     Ord e,
     Unbox e
   ) =>
@@ -934,16 +943,17 @@
   where
     Sz (m :. n) = size distMat
 
--- Make a search row for distances. Takes row x from a distance matrix and zips them with their
+-- | Make a search row for distances. Takes row x from a distance matrix and zips them with their
 -- column index. Then keeps only the valid elements of the row, that are still part of the available
 -- points. A minimum or maximum search can be performed on the resulting vector and a valid pair of
 -- distance and index can be obtained.
 searchRow ::
   ( PrimMonad m,
+    RealWorld ~ PrimState m,
     MonadThrow m,
     MonadUnliftIO m,
-    Mutable r Ix2 e,
-    Mutable r Ix1 (e, Ix1)
+    Mutable r e,
+    Mutable r (e, Ix1)
   ) =>
   Ix1 ->
   IntSet ->