diff --git a/demos/getting_started_1.py b/demos/getting_started_1.py
index b45e1de47a6e6f60453b0fd3fff75567d6a23085..8dbc52070e944786b8c78d8f56febf5e3f7ff767 100644
--- a/demos/getting_started_1.py
+++ b/demos/getting_started_1.py
@@ -93,12 +93,10 @@ if __name__ == '__main__':
 
     # Apply the mapping
     prior_correlation_structure = PD(ift.PS_field(power_space, power_spectrum))
-    print(prior_correlation_structure)
 
     # Insert the result into the diagonal of an harmonic space operator
     S = ift.DiagonalOperator(prior_correlation_structure)
     # S is the prior field covariance
-    assert False
 
     # Build instrument response consisting of a discretization, mask
     # and harmonic transformaion
diff --git a/demos/misc/convolution.py b/demos/misc/convolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..7d057775c54bd9acdce06da48fe15875876c735b
--- /dev/null
+++ b/demos/misc/convolution.py
@@ -0,0 +1,74 @@
+import numpy as np
+import nifty5 as ift
+
+
+def convtest(test_signal, delta, func):
+    domain = test_signal.domain
+
+    # Create Convolution Operator
+    conv_op = ift.FuncConvolutionOperator(domain, func)
+
+    # Convolve, Adjoint-Convolve
+    conv_signal = conv_op(test_signal)
+    cac_signal = conv_op.adjoint_times(conv_signal)
+
+    print(test_signal.integrate(), conv_signal.integrate(),
+          cac_signal.integrate())
+
+    # generate kernel image
+    conv_delta = conv_op(delta)
+
+    # Plot results
+    plot = ift.Plot()
+    plot.add(signal, title='Signal')
+    plot.add(conv_signal, title='Signal Convolved')
+    plot.add(cac_signal, title='Signal, Conv, Adj-Conv')
+    plot.add(conv_delta, title='Kernel')
+    plot.output()
+
+
+# Healpix test
+nside = 64
+npix = 12 * nside * nside
+
+domain = ift.HPSpace(nside)
+
+# Define test signal (some point sources)
+signal_vals = np.zeros(npix, dtype=np.float64)
+for i in range(0, npix, npix//12 + 27):
+    signal_vals[i] = 500.
+
+signal = ift.from_global_data(domain, signal_vals)
+
+delta_vals = np.zeros(npix, dtype=np.float64)
+delta_vals[0] = 1.0
+delta = ift.from_global_data(domain, delta_vals)
+
+
+# Define kernel function
+def func(theta):
+    ct = np.cos(theta)
+    return 1. * np.logical_and(ct > 0.7, ct <= 0.8)
+
+
+convtest(signal, delta, func)
+
+domain = ift.RGSpace((100, 100))
+# Define test signal (some point sources)
+signal_vals = np.zeros(domain.shape, dtype=np.float64)
+signal_vals[35, 70] = 5000.
+signal_vals[45, 8] = 5000.
+signal = ift.from_global_data(domain, signal_vals)
+
+# Define delta signal, generate kernel image
+delta_vals = np.zeros(domain.shape, dtype=np.float64)
+delta_vals[0, 0] = 1.0
+delta = ift.from_global_data(domain, delta_vals)
+
+
+# Define kernel function
+def func(dist):
+    return 1. * np.logical_and(dist > 0.1, dist <= 0.2)
+
+
+convtest(signal, delta, func)
diff --git a/docs/source/ift.rst b/docs/source/ift.rst
index 66acad68569d205d29946b3670a97060cca9bce7..3b297d1ee485c446ce1dff4b7c60aca45aa60647 100644
--- a/docs/source/ift.rst
+++ b/docs/source/ift.rst
@@ -37,7 +37,7 @@ NIFTy comes with reimplemented MAP and VI estimators.
 Free Theory & Implicit Operators
 --------------------------------
 
-A free IFT appears when the signal field :math:`{s}` and the noise :math:`{n}` of the data :math:`{d}` are independent, zero-centered Gaussian processes of kown covariances :math:`{S}` and :math:`{N}`, respectively,
+A free IFT appears when the signal field :math:`{s}` and the noise :math:`{n}` of the data :math:`{d}` are independent, zero-centered Gaussian processes of known covariances :math:`{S}` and :math:`{N}`, respectively,
 
 .. math::
 
@@ -49,7 +49,7 @@ and the measurement equation is linear in both signal and noise,
 
     d= R\, s + n,
 
-with :math:`{R}` being the measurement response, which maps the continous signal field into the discrete data space.
+with :math:`{R}` being the measurement response, which maps the continuous signal field into the discrete data space.
 
 This is called a free theory, as the information Hamiltonian
 
@@ -96,8 +96,8 @@ These implicit operators can be combined into new operators, e.g. to :math:`{D^{
 The invocation of an inverse operator applied to a vector might trigger the execution of a numerical linear algebra solver.
 
 Thus, when NIFTy calculates :math:`{m = D\, j}`, it actually solves :math:`{D^{-1} m = j}` for :math:`{m}` behind the scenes.
-The advantage of implicit operators to explicit matrices is the reduced memory requirements.
-The reconstruction of only a Megapixel image would otherwithe require the storage and processing of matrices with sizes of several Terabytes.
+The advantage of implicit operators compared to explicit matrices is the reduced memory consumption;
+for the reconstruction of just a Megapixel image the latter would already require several Terabytes.
 Larger images could not be dealt with due to the quadratic memory requirements of explicit operator representations.
 
 The demo codes `demos/getting_started_1.py` and `demos/Wiener_Filter.ipynb` illustrate this.
@@ -106,7 +106,7 @@ The demo codes `demos/getting_started_1.py` and `demos/Wiener_Filter.ipynb` illu
 Generative Models
 -----------------
 
-For more sophisticated measurement situations, involving non-linear measuremnts, unknown covariances, calibration constants and the like, it is recommended to formulate those as generative models for which NIFTy provides powerful inference algorithms.
+For more sophisticated measurement situations (involving non-linear measurements, unknown covariances, calibration constants and the like) it is recommended to formulate those as generative models for which NIFTy provides powerful inference algorithms.
 
 In a generative model, all known or unknown quantities are described as the results of generative processes, which start with simple probability distributions, like the uniform, the i.i.d. Gaussian, or the delta distribution.
 
@@ -129,7 +129,7 @@ NIFTy takes advantage of this formulation in several ways:
 1) All prior degrees of freedom have unit covariance, which improves the condition number of operators that need to be inverted.
 
 2) The amplitude operator can be regarded as part of the response, :math:`{R'=R\,A}`.
-   In general, more sophisticated responses can be constructed out of the composition of simpler operators.
+   In general, more sophisticated responses can be obtained by combining simpler operators.
 
 3) The response can be non-linear, e.g. :math:`{R'(s)=R \exp(A\,\xi)}`, see `demos/getting_started_2.py`.
 
@@ -168,7 +168,7 @@ Maximum a Posteriori
 
 One popular field estimation method is Maximum a Posteriori (MAP).
 
-It only requires to minimize the information Hamiltonian, e.g by a gradient descent method that stops when
+It only requires minimizing the information Hamiltonian, e.g. by a gradient descent method that stops when
 
 .. math::
 
diff --git a/nifty5/__init__.py b/nifty5/__init__.py
index a2581745e58ee702cb27b0e1be058f336e6449b4..c96769eccbb7411e2dd05e49f68326dbc51c22a3 100644
--- a/nifty5/__init__.py
+++ b/nifty5/__init__.py
@@ -50,6 +50,7 @@ from .operators.value_inserter import ValueInserter
 from .operators.energy_operators import (
     EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood,
     BernoulliEnergy, StandardHamiltonian, AveragedEnergy)
+from .operators.convolution_operators import FuncConvolutionOperator
 
 from .probing import probe_with_posterior_samples, probe_diagonal, \
     StatCalculator
@@ -81,7 +82,8 @@ from .library.dynamic_operator import (dynamic_operator,
 from .library.light_cone_operator import LightConeOperator
 
 from .library.wiener_filter_curvature import WienerFilterCurvature
-from .library.correlated_fields import CorrelatedField, MfCorrelatedField
+from .library.correlated_fields import (CorrelatedField, MfCorrelatedField,
+                                        MfPartiallyCorrelatedField)
 from .library.adjust_variances import (make_adjust_variances_hamiltonian,
                                        do_adjust_variances)
 
@@ -93,5 +95,10 @@ from .logger import logger
 
 from .linearization import Linearization
 
+from . import internal_config
+_scheme = internal_config.parallelization_scheme()
+if _scheme == "Samples":
+    from .minimization.metric_gaussian_kl_mpi import MetricGaussianKL_MPI
+
 # We deliberately don't set __all__ here, because we don't want people to do a
 # "from nifty5 import *"; that would swamp the global namespace.
diff --git a/nifty5/dobj.py b/nifty5/dobj.py
index 47dc229793b4a0c3aa37178f8effb6ec11081091..9f3d3dc7ea8b019fe9f3254b2d323dbd2d98a526 100644
--- a/nifty5/dobj.py
+++ b/nifty5/dobj.py
@@ -15,11 +15,19 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-try:
-    from mpi4py import MPI
-    if MPI.COMM_WORLD.Get_size() == 1:
-        from .data_objects.numpy_do import *
-    else:
-        from .data_objects.distributed_do import *
-except ImportError:
+from . import internal_config
+
+
+_scheme = internal_config.parallelization_scheme()
+
+if _scheme in ("Samples", "None"):
     from .data_objects.numpy_do import *
+else:
+    try:
+        from mpi4py import MPI
+        if MPI.COMM_WORLD.Get_size() == 1:
+            from .data_objects.numpy_do import *
+        else:
+            from .data_objects.distributed_do import *
+    except ImportError:
+        from .data_objects.numpy_do import *
diff --git a/nifty5/domains/lm_space.py b/nifty5/domains/lm_space.py
index f138e3e8803b09119cdb63f8efb038e6009f1377..5e0f30ffbba22267466895656855f9e4fed2de99 100644
--- a/nifty5/domains/lm_space.py
+++ b/nifty5/domains/lm_space.py
@@ -103,6 +103,36 @@ class LMSpace(StructuredDomain):
     def get_fft_smoothing_kernel_function(self, sigma):
         return lambda x: self._kernel(x, sigma)
 
+    def get_conv_kernel_from_func(self, func):
+        """Creates a convolution kernel defined by a function.
+
+        Parameters
+        ----------
+        func: function
+            This function needs to take exactly one argument, which is
+            colatitude in radians, and return the kernel amplitude at that
+            colatitude.
+
+        Assumes the function to be radially symmetric,
+        e.g. only dependant on theta in radians"""
+        from .gl_space import GLSpace
+        from ..operators.harmonic_operators import HarmonicTransformOperator
+        import pyHealpix
+        # define azimuthally symmetric spaces for kernel transform
+        gl = GLSpace(self.lmax + 1, 1)
+        lm0 = gl.get_default_codomain()
+        theta = pyHealpix.GL_thetas(gl.nlat)
+        # evaluate the kernel function at the required thetas
+        kernel_sphere = Field.from_global_data(gl, func(theta))
+        # normalize the kernel such that the integral over the sphere is 4pi
+        kernel_sphere = kernel_sphere * (4 * np.pi / kernel_sphere.integrate())
+        # compute the spherical harmonic coefficients of the kernel
+        op = HarmonicTransformOperator(lm0, gl)
+        kernel_lm = op.adjoint_times(kernel_sphere.weight(1)).to_global_data()
+        # evaluate the k lengths of the harmonic space
+        k_lengths = self.get_k_length_array().to_global_data().astype(np.int)
+        return Field.from_global_data(self, kernel_lm[k_lengths])
+
     @property
     def lmax(self):
         """int : maximum allowed :math:`l`
diff --git a/nifty5/domains/rg_space.py b/nifty5/domains/rg_space.py
index 2e1d2c0296c7a80bbdfe2d53bd9c849069323bee..02d4f106da9af7137632e8c738efe7bb8f4e9cc3 100644
--- a/nifty5/domains/rg_space.py
+++ b/nifty5/domains/rg_space.py
@@ -90,9 +90,7 @@ class RGSpace(StructuredDomain):
     def scalar_dvol(self):
         return self._dvol
 
-    def get_k_length_array(self):
-        if (not self.harmonic):
-            raise NotImplementedError
+    def _get_dist_array(self):
         ibegin = dobj.ibegin_from_shape(self._shape)
         res = np.arange(self.local_shape[0], dtype=np.float64) + ibegin[0]
         res = np.minimum(res, self.shape[0]-res)*self.distances[0]
@@ -106,6 +104,11 @@ class RGSpace(StructuredDomain):
             res = np.add.outer(res, tmp)
         return Field.from_local_data(self, np.sqrt(res))
 
+    def get_k_length_array(self):
+        if (not self.harmonic):
+            raise NotImplementedError
+        return self._get_dist_array()
+
     def get_unique_k_lengths(self):
         if (not self.harmonic):
             raise NotImplementedError
@@ -146,6 +149,27 @@ class RGSpace(StructuredDomain):
             raise NotImplementedError
         return lambda x: self._kernel(x, sigma)
 
+    def get_conv_kernel_from_func(self, func):
+        """Creates a convolution kernel defined by a function.
+
+        Parameters
+        ----------
+        func: function
+            This function needs to take exactly one argument, which is
+            distance from center (in the same units as the RGSpace distances),
+            and return the kernel amplitude at that distance.
+
+        Assumes the function to be radially symmetric,
+        e.g. only dependant on distance"""
+        from ..operators.harmonic_operators import HarmonicTransformOperator
+        if (not self.harmonic):
+            raise NotImplementedError
+        op = HarmonicTransformOperator(self, self.get_default_codomain())
+        dist = op.target[0]._get_dist_array()
+        kernel = Field.from_local_data(op.target, func(dist.local_data))
+        kernel = kernel / kernel.integrate()
+        return op.adjoint_times(kernel.weight(1))
+
     def get_default_codomain(self):
         """Returns a :class:`RGSpace` object representing the (position or
         harmonic) partner domain of `self`, depending on `self.harmonic`.
diff --git a/nifty5/field.py b/nifty5/field.py
index 79a8af7f581c02534f72c7b572154cb8bdd21b49..2929bcfa30ee7613282f0d4481314ab4c37d5b7b 100644
--- a/nifty5/field.py
+++ b/nifty5/field.py
@@ -636,6 +636,8 @@ class Field(object):
         return 0.5*(1.+self.tanh())
 
     def clip(self, min=None, max=None):
+        min = min.local_data if isinstance(min, Field) else min
+        max = max.local_data if isinstance(max, Field) else max
         return Field(self._domain, dobj.clip(self._val, min, max))
 
     def one_over(self):
diff --git a/nifty5/internal_config.py b/nifty5/internal_config.py
new file mode 100644
index 0000000000000000000000000000000000000000..c6c157aa0deb27bc40a2825aa99b04b30c893509
--- /dev/null
+++ b/nifty5/internal_config.py
@@ -0,0 +1,25 @@
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright(C) 2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+
+# Internal configuration switches, typically for experimental features.
+# Leave unchanged unless you know what you are doing!
+
+def parallelization_scheme():
+    return "Standard"
+    # return "Samples"
+    # return "None"
diff --git a/nifty5/library/correlated_fields.py b/nifty5/library/correlated_fields.py
index a349b15c48804937d068aedad6803799caef16a1..85297031d6fa1256fba83b4eadf12bdefe424ca1 100644
--- a/nifty5/library/correlated_fields.py
+++ b/nifty5/library/correlated_fields.py
@@ -141,3 +141,54 @@ def MfCorrelatedField(target, amplitudes, name='xi'):
     vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
 
     return ht(vol*A*ducktape(hsp, None, name))
+
+
+def MfPartiallyCorrelatedField(target, energy_amplitude, name='xi_u'):
+    """Constructs an operator which turns white Gaussian excitation fields
+    into a correlated field defined on a DomainTuple with two entries.
+    On the first domain, a white correlation structure is assumed. On
+    the second domain, the correlation structures given by energy_amplitude
+    is used.
+
+    This operator may be used as a model for multi-frequency reconstructions
+    with correlation structure only in the energy direction.
+
+    Parameters
+    ----------
+    target : Domain, DomainTuple or tuple of Domain
+        Target of the operator. Must contain exactly two spaces.
+        It is assumed that the second space is the energy domain.
+    energy_amplitude: Operator
+        amplitude operator for the energy correlation structure
+    name : string
+        :class:`MultiField` key for xi-field.
+
+    Returns
+    -------
+    Operator
+        Correlated field
+
+    Notes
+    -----
+    In NIFTy, non-harmonic RGSpaces are by definition periodic. Therefore
+    the operator constructed by this method will output a correlated field
+    with *periodic* boundary conditions. If a non-periodic field is needed,
+    one needs to combine this operator with a :class:`FieldZeroPadder` or even
+    two (one for the energy and one for the spatial subdomain)
+    """
+    tgt = DomainTuple.make(target)
+    if len(tgt) != 2:
+        raise ValueError
+
+    h_space = DomainTuple.make([dom.get_default_codomain() for dom in tgt])
+    ht1 = HarmonicTransformOperator(h_space, target=tgt[0], space=0)
+    ht2 = HarmonicTransformOperator(ht1.target, target=tgt[1], space=1)
+    ht = ht2 @ ht1
+
+    p_space = energy_amplitude.target[0]
+    power_distributor = PowerDistributor(h_space[-1], p_space)
+    A = power_distributor(energy_amplitude)
+
+    dd = ContractionOperator(h_space, 0).adjoint
+
+    return ht((dd @ A)*ducktape(h_space, None, name))
diff --git a/nifty5/library/smooth_linear_amplitude.py b/nifty5/library/smooth_linear_amplitude.py
index da5ad0c0af6e61e38819af3cf3238b0001e8c5cf..bdbe591ba49ee77a87d488acc81c6268c1795dcf 100644
--- a/nifty5/library/smooth_linear_amplitude.py
+++ b/nifty5/library/smooth_linear_amplitude.py
@@ -118,7 +118,7 @@ def CepstrumOperator(target, a, k0):
     return sym @ qht @ makeOp(cepstrum.sqrt())
 
 
-def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
+def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za=None, zq=None,
                 keys=['tau', 'phi', 'zero_mode']):
     '''Operator for parametrizing smooth amplitudes (square roots of power
     spectra).
@@ -168,10 +168,10 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
         LogRGSpace (see :class:`ExpTransform`).
     iv : float
         Prior standard deviation of y-intercept of power law.
-    za : float
+    za : float, optional
         The alpha-parameter of the inverse-gamma distribution.
         Setting the a seperate prior on the zeroGmode of the amplitude model.
-    zq : float
+    zq : float, optional
         The q-parameter of the inverse-gamma distribution.
 
     Returns
@@ -186,10 +186,17 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
 
     a, k0 = float(a), float(k0)
     sm, sv, im, iv = float(sm), float(sv), float(im), float(iv)
-    za, zq = float(za), float(zq)
     if sv <= 0 or iv <= 0:
         raise ValueError
 
+    if za is not None and zq is not None:
+        separate_zero_mode_prior = True
+        za, zq = float(za), float(zq)
+    else:
+        separate_zero_mode_prior = False
+        if za is not None or zq is not None:
+            raise ValueError("za and zq have to be given together")
+
     et = ExpTransform(target, n_pix)
     dom = et.domain[0]
 
@@ -208,14 +215,18 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
     # Combine linear and smooth component
     loglog_ampl = 0.5*(smooth + linear)
 
-    zero_mode = ValueInserter(et.target, (0,)*len(et.target.shape))
-
-    mask = np.ones(et.target.shape)
-    mask[(0,)*len(et.target.shape)] = 0.
-    mask = from_global_data(et.target, mask)
-    mask = DiagonalOperator(mask)
-    zero_mode = zero_mode @ InverseGammaOperator(
-        zero_mode.domain, za, zq) @ FieldAdapter(
-            zero_mode.domain, keys[2])
-    # Go from loglog-space to linear-linear-space
-    return (mask @ (et @ loglog_ampl.exp()) + zero_mode)
+    if not separate_zero_mode_prior:
+        # Go from loglog-space to linear-linear-space
+        return et @ loglog_ampl.exp()
+    else:
+        zero_mode = ValueInserter(et.target, (0,)*len(et.target.shape))
+
+        mask = np.ones(et.target.shape)
+        mask[(0,)*len(et.target.shape)] = 0.
+        mask = from_global_data(et.target, mask)
+        mask = DiagonalOperator(mask)
+        zero_mode = zero_mode @ InverseGammaOperator(
+            zero_mode.domain, za, zq) @ FieldAdapter(
+                zero_mode.domain, keys[2])
+
+        return mask @ (et @ loglog_ampl.exp()) + zero_mode
diff --git a/nifty5/minimization/metric_gaussian_kl_mpi.py b/nifty5/minimization/metric_gaussian_kl_mpi.py
new file mode 100644
index 0000000000000000000000000000000000000000..1970bb2e00cbd10858823acfce04c389a98ea632
--- /dev/null
+++ b/nifty5/minimization/metric_gaussian_kl_mpi.py
@@ -0,0 +1,149 @@
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright(C) 2013-2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+from .. import utilities
+from ..linearization import Linearization
+from ..operators.energy_operators import StandardHamiltonian
+from .energy import Energy
+from mpi4py import MPI
+import numpy as np
+from ..field import Field
+from ..multi_field import MultiField
+
+
+_comm = MPI.COMM_WORLD
+ntask = _comm.Get_size()
+rank = _comm.Get_rank()
+master = (rank == 0)
+
+
+def _shareRange(nwork, nshares, myshare):
+    nbase = nwork//nshares
+    additional = nwork % nshares
+    lo = myshare*nbase + min(myshare, additional)
+    hi = lo + nbase + int(myshare < additional)
+    return lo, hi
+
+
+def np_allreduce_sum(arr):
+    arr = np.array(arr)
+    res = np.empty_like(arr)
+    _comm.Allreduce(arr, res, MPI.SUM)
+    return res
+
+
+def allreduce_sum_field(fld):
+    if isinstance(fld, Field):
+        return Field.from_local_data(fld.domain,
+                                     np_allreduce_sum(fld.local_data))
+    res = tuple(
+        Field.from_local_data(f.domain, np_allreduce_sum(f.local_data))
+        for f in fld.values())
+    return MultiField(fld.domain, res)
+
+
+class MetricGaussianKL_MPI(Energy):
+    def __init__(self, mean, hamiltonian, n_samples, constants=[],
+                 point_estimates=[], mirror_samples=False,
+                 _samples=None):
+        super(MetricGaussianKL_MPI, self).__init__(mean)
+
+        if not isinstance(hamiltonian, StandardHamiltonian):
+            raise TypeError
+        if hamiltonian.domain is not mean.domain:
+            raise ValueError
+        if not isinstance(n_samples, int):
+            raise TypeError
+        self._constants = list(constants)
+        self._point_estimates = list(point_estimates)
+        if not isinstance(mirror_samples, bool):
+            raise TypeError
+
+        self._hamiltonian = hamiltonian
+
+        if _samples is None:
+            lo, hi = _shareRange(n_samples, ntask, rank)
+            met = hamiltonian(Linearization.make_partial_var(
+                mean, point_estimates, True)).metric
+            _samples = []
+            for i in range(lo, hi):
+                np.random.seed(i)
+                _samples.append(met.draw_sample(from_inverse=True))
+            if mirror_samples:
+                _samples += [-s for s in _samples]
+                n_samples *= 2
+            _samples = tuple(_samples)
+        self._samples = _samples
+        self._n_samples = n_samples
+        self._lin = Linearization.make_partial_var(mean, constants)
+        v, g = None, None
+        if len(self._samples) == 0:  # hack if there are too many MPI tasks
+            tmp = self._hamiltonian(self._lin)
+            v = 0. * tmp.val.local_data
+            g = 0. * tmp.gradient
+        else:
+            for s in self._samples:
+                tmp = self._hamiltonian(self._lin+s)
+                if v is None:
+                    v = tmp.val.local_data.copy()
+                    g = tmp.gradient
+                else:
+                    v += tmp.val.local_data
+                    g = g + tmp.gradient
+        self._val = np_allreduce_sum(v)[()] / self._n_samples
+        self._grad = allreduce_sum_field(g) / self._n_samples
+        self._metric = None
+
+    def at(self, position):
+        return MetricGaussianKL_MPI(
+            position, self._hamiltonian, self._n_samples, self._constants,
+            self._point_estimates, _samples=self._samples)
+
+    @property
+    def value(self):
+        return self._val
+
+    @property
+    def gradient(self):
+        return self._grad
+
+    def _get_metric(self):
+        lin = self._lin.with_want_metric()
+        if self._metric is None:
+            if len(self._samples) == 0:  # hack if there are too many MPI tasks
+                self._metric = self._hamiltonian(lin).metric.scale(0.)
+            else:
+                mymap = map(lambda v: self._hamiltonian(lin+v).metric,
+                            self._samples)
+                self._metric = utilities.my_sum(mymap)
+                self._metric = self._metric.scale(1./self._n_samples)
+
+    def apply_metric(self, x):
+        self._get_metric()
+        return allreduce_sum_field(self._metric(x))
+
+    @property
+    def metric(self):
+        if ntask > 1:
+            raise ValueError("not supported when MPI is active")
+        return self._metric
+
+    @property
+    def samples(self):
+        res = _comm.allgather(self._samples)
+        res = [item for sublist in res for item in sublist]
+        return res
diff --git a/nifty5/multi_field.py b/nifty5/multi_field.py
index 5206b85fb1049775072a9a8963c2a422a8c4c3a1..c9524f187a70465418ebc13ad3d043944a1586c3 100644
--- a/nifty5/multi_field.py
+++ b/nifty5/multi_field.py
@@ -192,8 +192,12 @@ class MultiField(object):
         return self._transform(lambda x: x.conjugate())
 
     def clip(self, min=None, max=None):
-        return MultiField(self._domain,
-                          tuple(clip(v, min, max) for v in self._val))
+        ncomp = len(self._val)
+        lmin = min._val if isinstance(min, MultiField) else (min,)*ncomp
+        lmax = max._val if isinstance(max, MultiField) else (max,)*ncomp
+        return MultiField(
+            self._domain,
+            tuple(self._val[i].clip(lmin[i], lmax[i]) for i in range(ncomp)))
 
     def all(self):
         for v in self._val:
diff --git a/nifty5/operators/convolution_operators.py b/nifty5/operators/convolution_operators.py
new file mode 100644
index 0000000000000000000000000000000000000000..26a2faea3956423bb229a935153169644c7212bd
--- /dev/null
+++ b/nifty5/operators/convolution_operators.py
@@ -0,0 +1,73 @@
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright(C) 2013-2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+import numpy as np
+
+from ..domains.rg_space import RGSpace
+from ..domains.lm_space import LMSpace
+from ..domains.hp_space import HPSpace
+from ..domains.gl_space import GLSpace
+from .endomorphic_operator import EndomorphicOperator
+from .harmonic_operators import HarmonicTransformOperator
+from .diagonal_operator import DiagonalOperator
+from .simple_linear_operators import WeightApplier
+from ..domain_tuple import DomainTuple
+from ..field import Field
+from .. import utilities
+
+
+def FuncConvolutionOperator(domain, func, space=None):
+    """Convolves input with a radially symmetric kernel defined by `func`
+
+    Parameters
+    ----------
+    domain: DomainTuple
+        Domain of the operator.
+    func: function
+        This function needs to take exactly one argument, which is
+        colatitude in radians, and return the kernel amplitude at that
+        colatitude.
+    space: int, optional
+        The index of the subdomain on which the operator should act
+        If None, it is set to 0 if `domain` contains exactly one space.
+        `domain[space]` must be of type `RGSpace`, `HPSpace`, or `GLSpace`.
+    """
+    domain = DomainTuple.make(domain)
+    space = utilities.infer_space(domain, space)
+    if not isinstance(domain[space], (RGSpace, HPSpace, GLSpace)):
+        raise TypeError("unsupported domain")
+    codomain = domain[space].get_default_codomain()
+    kernel = codomain.get_conv_kernel_from_func(func)
+    return _ConvolutionOperator(domain, kernel, space)
+
+
+def _ConvolutionOperator(domain, kernel, space=None):
+    domain = DomainTuple.make(domain)
+    space = utilities.infer_space(domain, space)
+    if len(kernel.domain) != 1:
+        raise ValueError("kernel needs exactly one domain")
+    if not isinstance(domain[space], (HPSpace, GLSpace, RGSpace)):
+        raise TypeError("need RGSpace, HPSpace, or GLSpace")
+    lm = [d for d in domain]
+    lm[space] = lm[space].get_default_codomain()
+    lm = DomainTuple.make(lm)
+    if lm[space] != kernel.domain[0]:
+        raise ValueError("Input domain and kernel are incompatible")
+    HT = HarmonicTransformOperator(lm, domain[space], space)
+    diag = DiagonalOperator(kernel*domain[space].total_volume, lm, (space,))
+    wgt = WeightApplier(domain, space, 1)
+    return HT(diag(HT.adjoint(wgt)))
diff --git a/nifty5/operators/simple_linear_operators.py b/nifty5/operators/simple_linear_operators.py
index 1628f2ab92cc15fe183d11023fb127e3c3bbdc47..3c1b54f4be5a62f7e7a51483a672b87746312d4f 100644
--- a/nifty5/operators/simple_linear_operators.py
+++ b/nifty5/operators/simple_linear_operators.py
@@ -63,6 +63,35 @@ class ConjugationOperator(EndomorphicOperator):
         return x.conjugate()
 
 
+class WeightApplier(EndomorphicOperator):
+    """Operator multiplying its input by a given power of dvol.
+
+    Parameters
+    ----------
+    domain: Domain, tuple of domains or DomainTuple
+        domain of the input field
+    spaces: list or tuple of int
+        indices of subdomains for which the weights shall be applied
+    power: int
+        the power of to be used for the volume factors
+
+    """
+    def __init__(self, domain, spaces, power):
+        from .. import utilities
+        self._domain = DomainTuple.make(domain)
+        if spaces is None:
+            self._spaces = None
+        else:
+            self._spaces = utilities.parse_spaces(spaces, len(self._domain))
+        self._power = int(power)
+        self._capability = self._all_ops
+
+    def apply(self, x, mode):
+        self._check_input(x, mode)
+        power = self._power if (mode & 3) else -self._power
+        return x.weight(power, spaces=self._spaces)
+
+
 class Realizer(EndomorphicOperator):
     """Operator returning the real component of its input.
 
diff --git a/test/test_operators/test_jacobian.py b/test/test_operators/test_jacobian.py
index 4bef7b25cd1e9c5ec202ad02214e925eee217fb8..cd0713251e96bdc40dfda3e04a62bd8186aa6f7b 100644
--- a/test/test_operators/test_jacobian.py
+++ b/test/test_operators/test_jacobian.py
@@ -87,6 +87,10 @@ def testBinary(type1, type2, space, seed):
     model = select_s1.clip(-1, 1)
     pos = ift.from_random("normal", dom1)
     ift.extra.check_jacobian_consistency(model, pos, ntries=20)
+    f = ift.from_random("normal", space)
+    model = select_s1.clip(f-0.1, f+1.)
+    pos = ift.from_random("normal", dom1)
+    ift.extra.check_jacobian_consistency(model, pos, ntries=20)
     if isinstance(space, ift.RGSpace):
         model = ift.FFTOperator(space)(select_s1*select_s2)
         pos = ift.from_random("normal", dom)