diff --git a/README.md b/README.md
index e846b3d..2efabe9 100644
--- a/README.md
+++ b/README.md
@@ -27,9 +27,9 @@ To install the latest stable version run
 
     pip install mpnum
 
-If you want to install `mpnum` from source, please run (on Unix)
+If you want to install this fork of `mpnum` from source, please run (on Unix)
 
-    git clone https://github.com/dseuss/mpnum.git
+    git clone https://github.com/moritzlange/mpnum.git
     cd mpnum
     pip install .
 
diff --git a/mpnum/mpsmpo.py b/mpnum/mpsmpo.py
index 337935d..8652229 100644
--- a/mpnum/mpsmpo.py
+++ b/mpnum/mpsmpo.py
@@ -125,16 +125,16 @@
 from __future__ import absolute_import, division, print_function
 
 import numpy as np
-from numpy.testing import assert_array_equal
 
 from six.moves import range
 
 from . import mparray as mp
 from .utils import local_to_global, matdot
 
+from scipy.linalg import eigh
 
 __all__ = ['mps_to_mpo', 'mps_to_pmps', 'pmps_dm_to_array',
-           'pmps_reduction', 'pmps_to_mpo', 'pmps_to_mps',
+           'pmps_reduction', 'mpo_to_pmps', 'pmps_to_mpo', 'pmps_to_mps',
            'reductions_mpo', 'reductions_mps_as_mpo',
            'reductions_mps_as_pmps', 'reductions_pmps', 'reductions']
 
@@ -243,7 +243,8 @@ def reductions_mpo(mpa, width=None, startsites=None, stopsites=None):
     startsites, stopsites = \
         _check_reductions_args(len(mpa), width, startsites, stopsites)
 
-    assert_array_equal(mpa.ndims, 2)
+    if not (np.array(mpa.ndims) == 2).all():
+        raise ValueError('Every site needs to have exactly two physical legs')
     rem_left = {0: np.array(1, ndmin=2)}
     rem_right = rem_left.copy()
 
@@ -372,6 +373,21 @@ def pmps_to_mpo(pmps):
     return mp.dot(pmps, pmps.adj())
 
 
+def mpo_to_pmps(mpo):
+    """Convert a tensor product MPO into a local purification MPS mixed state.
+
+    :param MPArray mpo: An MPA with two physical legs and rank one on all sites
+    :returns: An MPA with two physical legs (system and ancilla)
+
+    """
+    if not (np.array(mpo.ranks) == 1).all():
+        raise ValueError('Only implemented for rank-1 MPOs')
+    eigs = (eigh(lten[0, ..., 0] / np.trace(lten[0, ..., 0])) for lten in mpo.lt)
+    ltens = (np.sqrt(vals)[None, ...] * vecs for vals, vecs in eigs)
+    pmps = mp.MPArray.from_kron(ltens)
+    return pmps / mp.norm(pmps)
+
+
 def mps_to_pmps(mps):
     """Convert a pure MPS into a local purification MPS mixed state.
 
@@ -382,7 +398,8 @@ def mps_to_pmps(mps):
     :returns: An MPA with two physical legs (system and ancilla)
 
     """
-    assert_array_equal(mps.ndims, 1)
+    if not (np.array(mps.ndims) == 1).all():
+        raise ValueError('Only implemented for unit auxiliary dimensions')
     ltens = (lten.reshape(lten.shape[0:2] + (1, lten.shape[2])) for lten in mps.lt)
     return mp.MPArray(ltens)
 
diff --git a/tests/mpsmpo_test.py b/tests/mpsmpo_test.py
index 0e50b55..53b8858 100644
--- a/tests/mpsmpo_test.py
+++ b/tests/mpsmpo_test.py
@@ -5,6 +5,7 @@
 import numpy as np
 import pytest as pt
 from numpy.testing import assert_array_almost_equal
+from numpy.testing import assert_almost_equal
 
 import mpnum.factory as factory
 import mpnum.mparray as mp
@@ -223,6 +224,19 @@ def test_pmps_to_mpo(nr_sites, local_dim, rank, rgen):
     assert_array_almost_equal(rho_mp, rho_np)
 
 
+@pt.mark.parametrize('nr_sites, local_dim', [(1, 7), (2, 3), (3, 2), (6, 2),
+                                             (4, 3), (5, 2)])
+def test_mpo_to_pmps(nr_sites, local_dim, rgen):
+    pmps = factory.random_mpa(nr_sites, (local_dim, local_dim), 1,
+                              dtype=np.complex_, randstate=rgen)
+    assert_almost_equal(mp.normdist(mm.mpo_to_pmps(mm.pmps_to_mpo(pmps)),
+                                    pmps), 0)
+    mpo = factory.random_mpa(nr_sites, (local_dim, local_dim), 1,
+                             dtype=np.complex_, randstate=rgen, normalized=True)
+    assert_almost_equal(mp.normdist(mm.pmps_to_mpo(mm.mpo_to_pmps(mpo)),
+                                    mpo), 0)
+
+
 @pt.mark.parametrize('nr_sites, local_dim, rank', pt.MP_TEST_PARAMETERS)
 def test_mps_to_mpo(nr_sites, local_dim, rank, rgen):
     mps = factory.random_mps(nr_sites, local_dim, rank, randstate=rgen)