diff --git a/ChangeLog.md b/ChangeLog.md
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -5,9 +5,16 @@
 ## Unreleased changes
 
 
+## Version 0.7.0.1
+
+-   Random 1.2: Parallel functions now require an \`IOGenM\` random number
+    generator.
+-   Fix splitting of the random number generator.
+
+
 ## Version 0.6.1.1
 
--   Remove plotting functionality (gnuplot incompatible with ghc921).
+-   Remove plotting functionality (gnuplot incompatible with ghc922).
 -   Read files strictly.
 -   Refactor; flatten model hierarchy.
 
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
 
 # The ELynx Suite
 
-Version: 0.6.1.1.
+Version: 0.7.0.0.
 Reproducible evolution made easy.
 
 <p align="center"><img src="https://travis-ci.org/dschrempf/elynx.svg?branch=master"/></p>
@@ -74,9 +74,9 @@
     # OR: slynx --help
 
     Up to date
-    ELynx Suite version 0.6.1.1.
+    ELynx Suite version 0.7.0.0.
     Developed by Dominik Schrempf.
-    Compiled on February 22, 2022, at 15:10 pm, UTC.
+    Compiled on July 10, 2022, at 08:48 am, UTC.
     
     Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME]
                  [-f|--force] [--no-elynx-file] COMMAND
@@ -144,16 +144,16 @@
     # OR: slynx simulate --help
 
     Up to date
-    ELynx Suite version 0.6.1.1.
+    ELynx Suite version 0.7.0.0.
     Developed by Dominik Schrempf.
-    Compiled on February 22, 2022, at 15:10 pm, UTC.
+    Compiled on July 10, 2022, at 08:48 am, UTC.
     
     Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL]
                           [-m|--mixture-model MODEL] [-e|--edm-file NAME]
                           [-p|--siteprofile-files NAMES]
                           [-w|--mixture-model-weights "[DOUBLE,DOUBLE,...]"]
                           [-g|--gamma-rate-heterogeneity "(NCAT,SHAPE)"]
-                          (-l|--length NUMBER) [-S|--seed [INT]]
+                          (-l|--length NUMBER) [-S|--seed INT]
       Simulate multi sequence alignments.
     
     Available options:
@@ -175,8 +175,7 @@
       -g,--gamma-rate-heterogeneity "(NCAT,SHAPE)"
                                Number of gamma rate categories and shape parameter
       -l,--length NUMBER       Set alignment length to NUMBER
-      -S,--seed [INT]          Seed for random number generator; list of 32 bit
-                               integers with up to 256 elements (default: random)
+      -S,--seed INT            Seed for random number generator (default: random)
       -h,--help                Show this help text
     
     Substitution models:
diff --git a/bench/Bench.hs b/bench/Bench.hs
--- a/bench/Bench.hs
+++ b/bench/Bench.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  Main
 -- Description :  Benchmark Markov process along tree
--- Copyright   :  (c) 2021 Dominik Schrempf
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/elynx-markov.cabal b/elynx-markov.cabal
--- a/elynx-markov.cabal
+++ b/elynx-markov.cabal
@@ -1,55 +1,55 @@
-cabal-version:  3.0
-name:           elynx-markov
-version:        0.6.1.1
-synopsis:       Simulate molecular sequences along trees
-description:    Examine, modify, and simulate molecular sequences in a reproducible way. Please see the README on GitHub at <https://github.com/dschrempf/elynx>.
-category:       Bioinformatics
-homepage:       https://github.com/dschrempf/elynx#readme
-bug-reports:    https://github.com/dschrempf/elynx/issues
-author:         Dominik Schrempf
-maintainer:     dominik.schrempf@gmail.com
-copyright:      Dominik Schrempf (2021)
-license:        GPL-3.0-or-later
-license-file:   LICENSE
-build-type:     Simple
+cabal-version:      3.0
+name:               elynx-markov
+version:            0.7.0.1
+synopsis:           Simulate molecular sequences along trees
+description:
+  Examine, modify, and simulate molecular sequences in a reproducible way. Please see the README on GitHub at <https://github.com/dschrempf/elynx>.
 
+category:           Bioinformatics
+homepage:           https://github.com/dschrempf/elynx#readme
+bug-reports:        https://github.com/dschrempf/elynx/issues
+author:             Dominik Schrempf
+maintainer:         dominik.schrempf@gmail.com
+copyright:          2021 Dominik Schrempf
+license:            GPL-3.0-or-later
+license-file:       LICENSE
+build-type:         Simple
 extra-source-files:
-    README.md
-    ChangeLog.md
-    data/EDMDistsIQTREE.nex
-    data/EDMDistsPhylobayes.txt
-    data/HSSPMany.siteprofiles
-    data/HSSPSingle.siteprofiles
+  ChangeLog.md
+  data/EDMDistsIQTREE.nex
+  data/EDMDistsPhylobayes.txt
+  data/HSSPMany.siteprofiles
+  data/HSSPSingle.siteprofiles
+  README.md
 
 source-repository head
-  type: git
+  type:     git
   location: https://github.com/dschrempf/elynx
 
 library
   exposed-modules:
-      ELynx.MarkovProcess.AminoAcid
-      ELynx.MarkovProcess.CXXModels
-      ELynx.MarkovProcess.CXXModelsData
-      ELynx.MarkovProcess.GammaRateHeterogeneity
-      ELynx.MarkovProcess.MixtureModel
-      ELynx.MarkovProcess.Nucleotide
-      ELynx.MarkovProcess.PhyloModel
-      ELynx.MarkovProcess.RateMatrix
-      ELynx.MarkovProcess.SubstitutionModel
-      ELynx.Import.MarkovProcess.EDMModelPhylobayes
-      ELynx.Import.MarkovProcess.SiteprofilesPhylobayes
-      ELynx.Simulate.MarkovProcess
-      ELynx.Simulate.MarkovProcessAlongTree
-  other-modules:
-      Paths_elynx_markov
-  autogen-modules:
-      Paths_elynx_markov
-  hs-source-dirs: src
-  ghc-options: -Wall -Wunused-packages
+    ELynx.Import.MarkovProcess.EDMModelPhylobayes
+    ELynx.Import.MarkovProcess.SiteprofilesPhylobayes
+    ELynx.MarkovProcess.AminoAcid
+    ELynx.MarkovProcess.CXXModels
+    ELynx.MarkovProcess.CXXModelsData
+    ELynx.MarkovProcess.GammaRateHeterogeneity
+    ELynx.MarkovProcess.MixtureModel
+    ELynx.MarkovProcess.Nucleotide
+    ELynx.MarkovProcess.PhyloModel
+    ELynx.MarkovProcess.RateMatrix
+    ELynx.MarkovProcess.SubstitutionModel
+    ELynx.Simulate.MarkovProcess
+    ELynx.Simulate.MarkovProcessAlongTree
+
+  other-modules:    Paths_elynx_markov
+  autogen-modules:  Paths_elynx_markov
+  hs-source-dirs:   src
+  ghc-options:      -Wall -Wunused-packages
   build-depends:
-      async
+    , async
     , attoparsec
-    , base >=4.7 && <5
+    , base            >=4.7 && <5
     , bytestring
     , containers
     , elynx-seq
@@ -57,45 +57,47 @@
     , integration
     , math-functions
     , mwc-random
-    , primitive
+    , random
     , statistics
     , vector
+
   default-language: Haskell2010
 
 test-suite markov-test
-  type: exitcode-stdio-1.0
-  main-is: Spec.hs
+  type:             exitcode-stdio-1.0
+  main-is:          Spec.hs
   other-modules:
-      ELynx.MarkovProcess.AminoAcidSpec
-      ELynx.MarkovProcess.NucleotideSpec
-      ELynx.MarkovProcess.RateMatrixSpec
-      ELynx.Import.MarkovProcess.EDMModelPhylobayesSpec
-      ELynx.Import.MarkovProcess.SiteprofilesPhylobayesSpec
-      ELynx.Simulate.MarkovProcessAlongTreeSpec
-      Paths_elynx_markov
-  autogen-modules:
-      Paths_elynx_markov
-  hs-source-dirs: test
-  ghc-options: -Wall -Wunused-packages
+    ELynx.Import.MarkovProcess.EDMModelPhylobayesSpec
+    ELynx.Import.MarkovProcess.SiteprofilesPhylobayesSpec
+    ELynx.MarkovProcess.AminoAcidSpec
+    ELynx.MarkovProcess.NucleotideSpec
+    ELynx.MarkovProcess.RateMatrixSpec
+    ELynx.Simulate.MarkovProcessAlongTreeSpec
+    Paths_elynx_markov
+
+  autogen-modules:  Paths_elynx_markov
+  hs-source-dirs:   test
+  ghc-options:      -Wall -Wunused-packages
   build-depends:
-      base >=4.7 && <5
+    , base          >=4.7 && <5
     , containers
     , elynx-markov
     , elynx-tools
     , hmatrix
     , hspec
-    , mwc-random
+    , random
     , vector
+
   default-language: Haskell2010
 
 benchmark markov-bench
-  type:            exitcode-stdio-1.0
-  main-is:         Bench.hs
+  type:             exitcode-stdio-1.0
+  main-is:          Bench.hs
+
   -- other-modules:
   -- autogen-modules:
   -- ghc-options:
-  hs-source-dirs:  bench
-  ghc-options: -Wall -Wunused-packages
-  build-depends:
-      base
+  hs-source-dirs:   bench
+  ghc-options:      -Wall -Wunused-packages
+  build-depends:    base
   default-language: Haskell2010
diff --git a/src/ELynx/Import/MarkovProcess/EDMModelPhylobayes.hs b/src/ELynx/Import/MarkovProcess/EDMModelPhylobayes.hs
--- a/src/ELynx/Import/MarkovProcess/EDMModelPhylobayes.hs
+++ b/src/ELynx/Import/MarkovProcess/EDMModelPhylobayes.hs
@@ -3,7 +3,7 @@
 -- |
 -- Module      :  ELynx.Import.MarkovProcess.EDMModelPhylobayes
 -- Description :  Import stationary distributions from Phylobayes format
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs b/src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs
--- a/src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs
+++ b/src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  ELynx.Import.MarkovProcess.SiteprofilesPhylobayes
 -- Description :  Import site profiles in Phylobayes format
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/ELynx/MarkovProcess/AminoAcid.hs b/src/ELynx/MarkovProcess/AminoAcid.hs
--- a/src/ELynx/MarkovProcess/AminoAcid.hs
+++ b/src/ELynx/MarkovProcess/AminoAcid.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  ELynx.MarkovProcess.AminoAcid
 -- Description :  Amino acid rate matrices such as LG
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/ELynx/MarkovProcess/CXXModels.hs b/src/ELynx/MarkovProcess/CXXModels.hs
--- a/src/ELynx/MarkovProcess/CXXModels.hs
+++ b/src/ELynx/MarkovProcess/CXXModels.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  ELynx.MarkovProcess.CXXModels
 -- Description :  C10 to C60 models
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/ELynx/MarkovProcess/CXXModelsData.hs b/src/ELynx/MarkovProcess/CXXModelsData.hs
--- a/src/ELynx/MarkovProcess/CXXModelsData.hs
+++ b/src/ELynx/MarkovProcess/CXXModelsData.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  ELynx.MarkovProcess.CXXModelsData
 -- Description :  Stationary distributions and weights
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/ELynx/MarkovProcess/GammaRateHeterogeneity.hs b/src/ELynx/MarkovProcess/GammaRateHeterogeneity.hs
--- a/src/ELynx/MarkovProcess/GammaRateHeterogeneity.hs
+++ b/src/ELynx/MarkovProcess/GammaRateHeterogeneity.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  ELynx.MarkovProcess.GammaRateHeterogeneity
 -- Description :  Discrete gamma rate heterogeneity
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/ELynx/MarkovProcess/MixtureModel.hs b/src/ELynx/MarkovProcess/MixtureModel.hs
--- a/src/ELynx/MarkovProcess/MixtureModel.hs
+++ b/src/ELynx/MarkovProcess/MixtureModel.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  ELynx.MarkovProcess.MixtureModel
 -- Description :  Mixture models are a set of substitution models with weights
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/ELynx/MarkovProcess/Nucleotide.hs b/src/ELynx/MarkovProcess/Nucleotide.hs
--- a/src/ELynx/MarkovProcess/Nucleotide.hs
+++ b/src/ELynx/MarkovProcess/Nucleotide.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  ELynx.MarkovProcess.Nucleotide
 -- Description :  Substitution models using nucleotides
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/ELynx/MarkovProcess/PhyloModel.hs b/src/ELynx/MarkovProcess/PhyloModel.hs
--- a/src/ELynx/MarkovProcess/PhyloModel.hs
+++ b/src/ELynx/MarkovProcess/PhyloModel.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  ELynx.MarkovProcess.PhyloModel
 -- Description :  Phylogenetic model
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/ELynx/MarkovProcess/RateMatrix.hs b/src/ELynx/MarkovProcess/RateMatrix.hs
--- a/src/ELynx/MarkovProcess/RateMatrix.hs
+++ b/src/ELynx/MarkovProcess/RateMatrix.hs
@@ -2,7 +2,7 @@
 
 -- |
 -- Description :  Rate matrix helper functions
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPLv3
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -179,8 +179,8 @@
   | i == j = 0.0
   | i < j = es !! ijToKLower jI iI
   | otherwise =
-    error
-      "Float indices could not be compared during matrix creation."
+      error
+        "Float indices could not be compared during matrix creation."
   where
     iI = round i :: Int
     jI = round j :: Int
@@ -193,8 +193,8 @@
   | i == j = 0.0
   | i > j = es !! ijToKUpper n jI iI
   | otherwise =
-    error
-      "Float indices could not be compared during matrix creation."
+      error
+        "Float indices could not be compared during matrix creation."
   where
     iI = round i :: Int
     jI = round j :: Int
diff --git a/src/ELynx/MarkovProcess/SubstitutionModel.hs b/src/ELynx/MarkovProcess/SubstitutionModel.hs
--- a/src/ELynx/MarkovProcess/SubstitutionModel.hs
+++ b/src/ELynx/MarkovProcess/SubstitutionModel.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  ELynx.MarkovProcess.SubstitutionModel
 -- Description :  Data type describing substitution model
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/src/ELynx/Simulate/MarkovProcess.hs b/src/ELynx/Simulate/MarkovProcess.hs
--- a/src/ELynx/Simulate/MarkovProcess.hs
+++ b/src/ELynx/Simulate/MarkovProcess.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  ELynx.Simulate.MarkovProcess
 -- Description :  Markov process helpers
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -17,11 +17,10 @@
   )
 where
 
-import Control.Monad.Primitive
 import ELynx.MarkovProcess.RateMatrix
 import Numeric.LinearAlgebra
-import System.Random.MWC
 import System.Random.MWC.Distributions
+import System.Random.Stateful
 
 -- | A probability matrix, P_ij(t) = Pr (X_t = j | X_0 = i).
 type ProbMatrix = Matrix R
@@ -33,19 +32,20 @@
 -- another in a specific time (branch length).
 probMatrix :: RateMatrix -> Double -> ProbMatrix
 probMatrix q t
-  | t == 0 =
-    if rows q == cols q
-      then ident (rows q)
-      else error "probMatrix: Matrix is not square."
+  | n /= m = error "probMatrix: Matrix is not square."
+  | t == 0 = ident n
   | t < 0 = error "probMatrix: Time is negative."
   | otherwise = expm $ scale t q
+  where
+    n = rows q
+    m = cols q
 
 -- | Move from a given state to a new one according to a transition probability
 -- matrix .
 --
 -- This function is the bottleneck of the simulator and takes up most of the
 -- computation time.
-jump :: (PrimMonad m) => State -> ProbMatrix -> Gen (PrimState m) -> m State
+jump :: StatefulGen g m => State -> ProbMatrix -> g -> m State
 jump i p = categorical (p ! i)
 
 -- XXX: Maybe for later, use condensed tables.
@@ -55,7 +55,7 @@
 --
 -- Do not generate table for each jump.
 --
--- jump :: (PrimMonad m) => State -> ProbMatrix -> Gen (PrimState m) -> m State
+-- jump :: StatefulGen g m => State -> ProbMatrix -> g -> m State
 -- jump i p = genFromTable table
 --   where
 --     ws = toList $ p ! i
diff --git a/src/ELynx/Simulate/MarkovProcessAlongTree.hs b/src/ELynx/Simulate/MarkovProcessAlongTree.hs
--- a/src/ELynx/Simulate/MarkovProcessAlongTree.hs
+++ b/src/ELynx/Simulate/MarkovProcessAlongTree.hs
@@ -2,7 +2,7 @@
 
 -- |
 --   Description :  Work with transition probability matrices on rooted trees
---   Copyright   :  (c) Dominik Schrempf 2021
+--   Copyright   :  2021 Dominik Schrempf
 --   License     :  GPLv3
 --
 --   Maintainer  :  dominik.schrempf@gmail.com
@@ -30,14 +30,12 @@
 import Control.Concurrent
 import Control.Concurrent.Async
 import Control.Monad
-import Control.Monad.Primitive
 import Data.Tree
 import qualified Data.Vector as V
-import Data.Word (Word32)
 import ELynx.MarkovProcess.RateMatrix
 import ELynx.Simulate.MarkovProcess
-import System.Random.MWC
-import System.Random.MWC.Distributions (categorical)
+import System.Random.MWC.Distributions
+import System.Random.Stateful
 
 -- XXX @performace. The horizontal concatenation might be slow. If so,
 -- 'concatenateSeqs' or 'concatenateAlignments' can be used, which directly
@@ -55,10 +53,10 @@
 toProbTree q = fmap (probMatrix q)
 
 getRootStates ::
-  PrimMonad m =>
+  StatefulGen g m =>
   Int ->
   StationaryDistribution ->
-  Gen (PrimState m) ->
+  g ->
   m [State]
 getRootStates n d g = replicateM n $ categorical d g
 
@@ -69,12 +67,12 @@
 -- XXX: Improve performance. Use vectors, not lists. I am actually not sure if
 -- this improves performance...
 simulateAndFlatten ::
-  PrimMonad m =>
+  StatefulGen g m =>
   Int ->
   StationaryDistribution ->
   ExchangeabilityMatrix ->
   Tree Double ->
-  Gen (PrimState m) ->
+  g ->
   m [[State]]
 simulateAndFlatten n d e t g = do
   let q = fromExchangeabilityMatrix e d
@@ -86,10 +84,10 @@
 -- Recursively jump down the branches to the leafs. Forget states at internal
 -- nodes.
 simulateAndFlatten' ::
-  (PrimMonad m) =>
+  StatefulGen g m =>
   [State] ->
   Tree ProbMatrix ->
-  Gen (PrimState m) ->
+  g ->
   m [[State]]
 simulateAndFlatten' is (Node p f) g = do
   is' <- mapM (\i -> jump i p g) is
@@ -99,15 +97,17 @@
 
 -- | See 'simulateAndFlatten', parallel version.
 simulateAndFlattenPar ::
+  RandomGen g =>
   Int ->
   StationaryDistribution ->
   ExchangeabilityMatrix ->
   Tree Double ->
-  GenIO ->
+  IOGenM g ->
   IO [[State]]
 simulateAndFlattenPar n d e t g = do
   c <- getNumCapabilities
-  gs <- splitGen c g
+  rs <- replicateM c $ splitGenM g
+  gs <- mapM newIOGenM rs
   let chunks = getChunks c n
       q = fromExchangeabilityMatrix e d
       pt = toProbTree q t
@@ -126,12 +126,12 @@
 -- internal nodes. The result is a tree with the list of simulated states as
 -- node labels.
 simulate ::
-  PrimMonad m =>
+  StatefulGen g m =>
   Int ->
   StationaryDistribution ->
   ExchangeabilityMatrix ->
   Tree Double ->
-  Gen (PrimState m) ->
+  g ->
   m (Tree [State])
 simulate n d e t g = do
   let q = fromExchangeabilityMatrix e d
@@ -142,10 +142,10 @@
 -- This is the heart of the simulation. Take a tree and a list of root states.
 -- Recursively jump down the branches to the leafs.
 simulate' ::
-  (PrimMonad m) =>
+  StatefulGen g m =>
   [State] ->
   Tree ProbMatrix ->
-  Gen (PrimState m) ->
+  g ->
   m (Tree [State])
 simulate' is (Node p f) g = do
   is' <- mapM (\i -> jump i p g) is
@@ -159,11 +159,11 @@
   fmap (\a -> V.map (`probMatrix` a) qs)
 
 getComponentsAndRootStates ::
-  PrimMonad m =>
+  StatefulGen g m =>
   Int ->
   V.Vector Double ->
   V.Vector StationaryDistribution ->
-  Gen (PrimState m) ->
+  g ->
   m ([Int], [State])
 getComponentsAndRootStates n ws ds g = do
   cs <- replicateM n $ categorical ws g
@@ -174,13 +174,13 @@
 -- corresponding weights. Forget states at internal nodes. See also
 -- 'simulateAndFlatten'.
 simulateAndFlattenMixtureModel ::
-  PrimMonad m =>
+  StatefulGen g m =>
   Int ->
   V.Vector Double ->
   V.Vector StationaryDistribution ->
   V.Vector ExchangeabilityMatrix ->
   Tree Double ->
-  Gen (PrimState m) ->
+  g ->
   -- | (IndicesOfComponents, [SimulatedSequenceForEachTip])
   m ([Int], [[State]])
 simulateAndFlattenMixtureModel n ws ds es t g = do
@@ -191,11 +191,11 @@
   return (cs, ss)
 
 simulateAndFlattenMixtureModel' ::
-  (PrimMonad m) =>
+  StatefulGen g m =>
   [State] ->
   [Int] ->
   Tree (V.Vector ProbMatrix) ->
-  Gen (PrimState m) ->
+  g ->
   m [[State]]
 simulateAndFlattenMixtureModel' is cs (Node ps f) g = do
   is' <- sequence [jump i (ps V.! c) g | (i, c) <- zip is cs]
@@ -212,28 +212,27 @@
     r = n `mod` c
     ns = replicate r (n' + 1) ++ replicate (c - r) n'
 
-splitGen :: PrimMonad m => Int -> Gen (PrimState m) -> m [Gen (PrimState m)]
-splitGen n gen
-  | n <= 0 = return []
-  | otherwise = do
-    seeds :: [V.Vector Word32] <- replicateM (n -1) $ uniformVector gen 256
-    fmap (gen :) (mapM initialize seeds)
-
-parComp :: Int -> (Int -> GenIO -> IO b) -> GenIO -> IO [b]
+-- NOTE: We could move away from IO here, but moving away from 'mapConcurrently'
+-- requires benchmarks. I am just not sure if it makes sense to spend more time
+-- on this since the parallelization itself is a bit weird. Like so, we walk
+-- along separate trees in each process.
+parComp :: RandomGen g => Int -> (Int -> IOGenM g -> IO b) -> IOGenM g -> IO [b]
 parComp num fun gen = do
   ncap <- getNumCapabilities
   let chunks = getChunks ncap num
-  gs <- splitGen ncap gen
+  rs <- replicateM ncap $ splitGenM gen
+  gs <- mapM newIOGenM rs
   mapConcurrently (uncurry fun) (zip chunks gs)
 
 -- | See 'simulateAndFlattenMixtureModel', parallel version.
 simulateAndFlattenMixtureModelPar ::
+  RandomGen g =>
   Int ->
   V.Vector Double ->
   V.Vector StationaryDistribution ->
   V.Vector ExchangeabilityMatrix ->
   Tree Double ->
-  GenIO ->
+  IOGenM g ->
   IO ([Int], [[State]])
 simulateAndFlattenMixtureModelPar n ws ds es t g = do
   let qs = V.zipWith fromExchangeabilityMatrix es ds
@@ -256,13 +255,13 @@
 -- corresponding weights. Keep states at internal nodes. See also
 -- 'simulate'.
 simulateMixtureModel ::
-  PrimMonad m =>
+  StatefulGen g m =>
   Int ->
   V.Vector Double ->
   V.Vector StationaryDistribution ->
   V.Vector ExchangeabilityMatrix ->
   Tree Double ->
-  Gen (PrimState m) ->
+  g ->
   m (Tree [State])
 simulateMixtureModel n ws ds es t g = do
   let qs = V.zipWith fromExchangeabilityMatrix es ds
@@ -273,11 +272,11 @@
 -- See 'simulateAlongProbTree', only we have a number of mixture components. The
 -- starting states and the components for each site have to be provided.
 simulateMixtureModel' ::
-  (PrimMonad m) =>
+  StatefulGen g m =>
   [State] ->
   [Int] ->
   Tree (V.Vector ProbMatrix) ->
-  Gen (PrimState m) ->
+  g ->
   m (Tree [State])
 simulateMixtureModel' is cs (Node ps f) g = do
   is' <- sequence [jump i (ps V.! c) g | (i, c) <- zip is cs]
diff --git a/test/ELynx/Import/MarkovProcess/EDMModelPhylobayesSpec.hs b/test/ELynx/Import/MarkovProcess/EDMModelPhylobayesSpec.hs
--- a/test/ELynx/Import/MarkovProcess/EDMModelPhylobayesSpec.hs
+++ b/test/ELynx/Import/MarkovProcess/EDMModelPhylobayesSpec.hs
@@ -1,6 +1,6 @@
 -- |
 -- Module      :  ELynx.Import.MarkovProcess.EDMModelPhylobayesSpec
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/test/ELynx/Import/MarkovProcess/SiteprofilesPhylobayesSpec.hs b/test/ELynx/Import/MarkovProcess/SiteprofilesPhylobayesSpec.hs
--- a/test/ELynx/Import/MarkovProcess/SiteprofilesPhylobayesSpec.hs
+++ b/test/ELynx/Import/MarkovProcess/SiteprofilesPhylobayesSpec.hs
@@ -1,6 +1,6 @@
 -- |
 -- Module      :  ELynx.Import.MarkovProcess.SiteprofilesPhylobayesSpec
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/test/ELynx/MarkovProcess/AminoAcidSpec.hs b/test/ELynx/MarkovProcess/AminoAcidSpec.hs
--- a/test/ELynx/MarkovProcess/AminoAcidSpec.hs
+++ b/test/ELynx/MarkovProcess/AminoAcidSpec.hs
@@ -1,6 +1,6 @@
 -- |
 -- Module      :  ELynx.MarkovProcess.AminoAcidSpec
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/test/ELynx/MarkovProcess/NucleotideSpec.hs b/test/ELynx/MarkovProcess/NucleotideSpec.hs
--- a/test/ELynx/MarkovProcess/NucleotideSpec.hs
+++ b/test/ELynx/MarkovProcess/NucleotideSpec.hs
@@ -1,6 +1,6 @@
 -- |
 -- Module      :  ELynx.MarkovProcess.NucleotideSpec
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/test/ELynx/MarkovProcess/RateMatrixSpec.hs b/test/ELynx/MarkovProcess/RateMatrixSpec.hs
--- a/test/ELynx/MarkovProcess/RateMatrixSpec.hs
+++ b/test/ELynx/MarkovProcess/RateMatrixSpec.hs
@@ -1,7 +1,7 @@
 -- |
 -- Module      :  ELynx.MarkovProcess.RateMatrixSpec
 -- Description :  Unit tests for rate matrices
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
diff --git a/test/ELynx/Simulate/MarkovProcessAlongTreeSpec.hs b/test/ELynx/Simulate/MarkovProcessAlongTreeSpec.hs
--- a/test/ELynx/Simulate/MarkovProcessAlongTreeSpec.hs
+++ b/test/ELynx/Simulate/MarkovProcessAlongTreeSpec.hs
@@ -1,6 +1,6 @@
 -- |
 -- Module      :  ELynx.Simulate.MarkovProcessAlongTreeSpec
--- Copyright   :  (c) Dominik Schrempf 2021
+-- Copyright   :  2021 Dominik Schrempf
 -- License     :  GPL-3.0-or-later
 --
 -- Maintainer  :  dominik.schrempf@gmail.com
@@ -18,7 +18,7 @@
 import qualified ELynx.MarkovProcess.SubstitutionModel as S
 import ELynx.Simulate.MarkovProcess
 import ELynx.Simulate.MarkovProcessAlongTree
-import System.Random.MWC
+import System.Random.Stateful
 import Test.Hspec
 
 -- sampleNewickText :: T.Text
@@ -42,11 +42,11 @@
 spec :: Spec
 spec = describe "simulateNSitesAlongTree" $ do
   it "simulates one site along an easy tree" $ do
-    gen <- create
+    gen <- newIOGenM $ mkStdGen 0
     tr <- simulate 1 d e nullTree gen
     fmap head tr `shouldBe` nullStateTree
   it "simulates some sites along a harder tree" $ do
-    gen <- create
+    gen <- newIOGenM $ mkStdGen 0
     tr <- simulate 10 d e testTree gen
     (length . rootLabel $ tr) `shouldBe` 10
   where