diff --git a/demos/getting_started_1.py b/demos/getting_started_1.py
index 4fe211864d9f60edb04e4cfe192ea2b088ff5cd1..54399816fce70cab4b2adf12df390cd3e9fc543e 100644
--- a/demos/getting_started_1.py
+++ b/demos/getting_started_1.py
@@ -6,19 +6,19 @@ def make_chess_mask():
     mask = np.ones(position_space.shape)
     for i in range(4):
         for j in range(4):
-            if (i+j)%2 == 0:
+            if (i+j) % 2 == 0:
                 mask[i*128//4:(i+1)*128//4, j*128//4:(j+1)*128//4] = 0
     return mask
 
 
 def make_random_mask():
-    mask = ift.from_random('pm1',position_space)
+    mask = ift.from_random('pm1', position_space)
     mask = (mask+1)/2
     return mask.val
 
 
 if __name__ == '__main__':
-    ## describtion of the tutorial ###
+    # # description of the tutorial ###
 
     # Choose problem geometry and masking
 
@@ -34,11 +34,11 @@ if __name__ == '__main__':
     # position_space = ift.HPSpace(128)
     # mask = make_random_mask()
 
-
     harmonic_space = position_space.get_default_codomain()
     HT = ift.HarmonicTransformOperator(harmonic_space, target=position_space)
 
-    # set correlation structure with a power spectrum and build prior correlation covariance
+    # set correlation structure with a power spectrum and build
+    # prior correlation covariance
     def power_spectrum(k):
         return 100. / (20.+k**3)
     power_space = ift.PowerSpace(harmonic_space)
@@ -47,9 +47,10 @@ if __name__ == '__main__':
 
     S = ift.DiagonalOperator(prior_correlation_structure)
 
-    # build instrument response consisting of a discretization, mask and harmonic transformaion
+    # build instrument response consisting of a discretization, mask
+    # and harmonic transformaion
     GR = ift.GeometryRemover(position_space)
-    mask = ift.Field(position_space,val=mask)
+    mask = ift.Field.from_global_data(position_space, mask)
     Mask = ift.DiagonalOperator(mask)
     R = GR * Mask * HT
 
@@ -69,10 +70,10 @@ if __name__ == '__main__':
     D_inv = R.adjoint * N.inverse * R + S.inverse
     # make it invertible
     IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3)
-    D = ift.InversionEnabler(D_inv,IC,approximation=S.inverse).inverse
+    D = ift.InversionEnabler(D_inv, IC, approximation=S.inverse).inverse
 
     # WIENER FILTER
     m = D(j)
 
-    ##PLOTTING
-    #Truth, data, reconstruction, residuals
+    # PLOTTING
+    # Truth, data, reconstruction, residuals
diff --git a/demos/getting_started_2.py b/demos/getting_started_2.py
index c0f2b773adfab51903d18f643931101ca1e3c630..0fa75285be5ab744056a74761b1eb8d8b15d76e3 100644
--- a/demos/getting_started_2.py
+++ b/demos/getting_started_2.py
@@ -6,14 +6,14 @@ def get_2D_exposure():
     x_shape, y_shape = position_space.shape
 
     exposure = np.ones(position_space.shape)
-    exposure[x_shape/3:x_shape/2,:] *= 2.
-    exposure[x_shape*4/5:x_shape,:] *= .1
-    exposure[x_shape/2:x_shape*3/2,:] *=3.
-    exposure[:,x_shape / 3:x_shape / 2] *= 2.
-    exposure[:,x_shape * 4 / 5:x_shape] *= .1
-    exposure[:,x_shape / 2:x_shape * 3 / 2] *= 3.
-
-    exposure = ift.Field(position_space, val=exposure)
+    exposure[x_shape/3:x_shape/2, :] *= 2.
+    exposure[x_shape*4/5:x_shape, :] *= .1
+    exposure[x_shape/2:x_shape*3/2, :] *= 3.
+    exposure[:, x_shape/3:x_shape/2] *= 2.
+    exposure[:, x_shape*4/5:x_shape] *= .1
+    exposure[:, x_shape/2:x_shape*3/2] *= 3.
+
+    exposure = ift.Field.from_global_data(position_space, exposure)
     return exposure
 
 
@@ -27,14 +27,13 @@ if __name__ == '__main__':
     # position_space = ift.RGSpace(1024)
     # exposure = np.ones(position_space.shape)
 
-
-    # Two dimensional regular grid with inhomogeneous exposure
+    # Two-dimensional regular grid with inhomogeneous exposure
     position_space = ift.RGSpace([512, 512])
     exposure = get_2D_exposure()
 
     # # Sphere with with uniform exposure
     # position_space = ift.HPSpace(128)
-    # exposure = np.ones(position_space.shape)
+    # exposure = ift.Field.full(position_space, 1.)
 
     # Defining harmonic space and transform
     harmonic_space = position_space.get_default_codomain()
@@ -58,7 +57,6 @@ if __name__ == '__main__':
     logsky = HT(logsky_h)
     sky = ift.PointwiseExponential(logsky)
 
-    exposure = ift.Field(position_space, val=exposure)
     M = ift.DiagonalOperator(exposure)
     GR = ift.GeometryRemover(position_space)
     # Set up instrumental response
@@ -68,14 +66,16 @@ if __name__ == '__main__':
     d_space = R.target[0]
     lamb = R(sky)
     mock_position = ift.from_random('normal', lamb.position.domain)
-    data = np.random.poisson(lamb.at(mock_position).value.val.astype(np.float64))
-    data = ift.Field.from_local_data(d_space, data)
+    data = lamb.at(mock_position).value
+    data = np.random.poisson(data.to_global_data().astype(np.float64))
+    data = ift.Field.from_global_data(d_space, data)
 
     # Compute likelihood and Hamiltonian
     position = ift.from_random('normal', lamb.position.domain)
     likelihood = ift.PoissonianEnergy(lamb, data)
     ic_cg = ift.GradientNormController(iteration_limit=50)
-    ic_newton = ift.GradientNormController(name='Newton',iteration_limit=50, tol_abs_gradnorm=1e-3)
+    ic_newton = ift.GradientNormController(name='Newton', iteration_limit=50,
+                                           tol_abs_gradnorm=1e-3)
     minimizer = ift.RelaxedNewton(ic_newton)
 
     # Minimize the Hamiltonian
@@ -84,6 +84,4 @@ if __name__ == '__main__':
 
     # Plot results
     result_sky = sky.at(H.position).value
-    ##PLOTTING
-
-
+    # PLOTTING
diff --git a/demos/getting_started_3.py b/demos/getting_started_3.py
index e9c670e9fa0ec87c72feeb4367f529e2738e9885..1952cd0a560eb6b932529632485a5ed881c72e8b 100644
--- a/demos/getting_started_3.py
+++ b/demos/getting_started_3.py
@@ -7,18 +7,20 @@ from scipy.io import loadmat
 
 
 def get_random_LOS(n_los):
-    starts = list(np.random.uniform(0,1,(n_los,2)).T)
-    ends = list(np.random.uniform(0,1,(n_los,2)).T)
+    starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
+    ends = list(np.random.uniform(0, 1, (n_los, 2)).T)
 
     return starts, ends
 
+
 if __name__ == '__main__':
-    ### ABOUT THIS TUTORIAL
+    # ## ABOUT THIS TUTORIAL
     np.random.seed(42)
-    position_space = ift.RGSpace([128,128])
+    position_space = ift.RGSpace([128, 128])
 
     # Setting up an amplitude model
-    A, amplitude_internals = make_amplitude_model(position_space,16, 1, 10, -4., 1, 0., 1.)
+    A, amplitude_internals = make_amplitude_model(
+        position_space, 16, 1, 10, -4., 1, 0., 1.)
 
     # Building the model for a correlated signal
     harmonic_space = position_space.get_default_codomain()
@@ -46,8 +48,7 @@ if __name__ == '__main__':
     # specify noise
     data_space = R.target
     noise = .001
-    N = ift.ScalingOperator(noise,data_space)
-
+    N = ift.ScalingOperator(noise, data_space)
 
     # generate mock data
     MOCK_POSITION = ift.from_random('normal', signal.position.domain)
@@ -63,28 +64,31 @@ if __name__ == '__main__':
     minimizer = ift.RelaxedNewton(ic_newton)
 
     # build model Hamiltonian
-    H = ift.Hamiltonian(likelihood,ic_cg,iteration_controller_sampling=ic_sampling)
+    H = ift.Hamiltonian(likelihood, ic_cg,
+                        iteration_controller_sampling=ic_sampling)
 
-    INITIAL_POSITION = ift.from_random('normal',H.position.domain)
+    INITIAL_POSITION = ift.from_random('normal', H.position.domain)
     position = INITIAL_POSITION
 
-    ift.plot(signal.at(MOCK_POSITION).value,name='truth.pdf')
-    ift.plot(R.adjoint_times(data),name='data.pdf')
-    ift.plot([ A.at(MOCK_POSITION).value], name='power.pdf')
+    ift.plot(signal.at(MOCK_POSITION).value, name='truth.pdf')
+    ift.plot(R.adjoint_times(data), name='data.pdf')
+    ift.plot([A.at(MOCK_POSITION).value], name='power.pdf')
 
     # number of samples used to estimate the KL
     N_samples = 20
     for i in range(5):
         H = H.at(position)
-        samples = [H.curvature.draw_sample(from_inverse=True) for _ in range(N_samples)]
+        samples = [H.curvature.draw_sample(from_inverse=True)
+                   for _ in range(N_samples)]
 
         KL = ift.SampledKullbachLeiblerDivergence(H, samples, ic_cg)
         KL, convergence = minimizer(KL)
         position = KL.position
 
-        ift.plot(signal.at(position).value,name='reconstruction.pdf')
+        ift.plot(signal.at(position).value, name='reconstruction.pdf')
 
-        ift.plot([A.at(position).value, A.at(MOCK_POSITION).value],name='power.pdf')
+        ift.plot([A.at(position).value, A.at(MOCK_POSITION).value],
+                 name='power.pdf')
 
     avrg = 0.
     va = 0.
@@ -101,15 +105,5 @@ if __name__ == '__main__':
     std = ift.sqrt(va)
     ift.plot(avrg, name='avrg.pdf')
     ift.plot(std, name='std.pdf')
-    ift.plot([A.at(position).value, A.at(MOCK_POSITION).value]+powers, name='power.pdf')
-
-
-
-
-
-
-
-
-
-
-
+    ift.plot([A.at(position).value, A.at(MOCK_POSITION).value]+powers,
+             name='power.pdf')
diff --git a/nifty5/library/amplitude_model.py b/nifty5/library/amplitude_model.py
index a85eace4292cd3dc9bc4e08045dc86148846352b..630b1678b17e9de1c10420074e7b2c047bf46d58 100644
--- a/nifty5/library/amplitude_model.py
+++ b/nifty5/library/amplitude_model.py
@@ -22,7 +22,7 @@ def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
 
     Npixdof : #pix in dof_space
 
-    ceps_a, ceps_k0 : Smoothnessparameters in ceps_kernel
+    ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel
                         eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
                         a = ceps_a,  k0 = ceps_k0
 
@@ -48,7 +48,7 @@ def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
     phi_sig = np.array([sv, iv])
 
     slope = SlopeOperator(param_space, logk_space, phi_sig)
-    norm_phi_mean = Field(param_space, val=phi_mean/phi_sig)
+    norm_phi_mean = Field.from_global_data(param_space, phi_mean/phi_sig)
 
     fields = {keys[0]: Field.from_random('normal', dof_space),
               keys[1]: Field.from_random('normal', param_space)}
@@ -93,7 +93,7 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
     # Prepare q_array
     q_array = np.zeros((dim,) + shape)
     if dim == 1:
-        ks = domain.get_k_length_array().val
+        ks = domain.get_k_length_array().to_global_data()
         q_array = np.array([ks])
     else:
         for i in range(dim):
@@ -120,4 +120,4 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
         # Do summation
         cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
 
-    return Field(domain, val=cepstrum_field)
+    return Field.from_global_data(domain, cepstrum_field)
diff --git a/nifty5/multi/multi_field.py b/nifty5/multi/multi_field.py
index f2257337a19c4c628c9a3e9bdd3b5c6225954d8a..f8074b03c15eaa75b4ae75b8ae0d015bce4a4cce 100644
--- a/nifty5/multi/multi_field.py
+++ b/nifty5/multi/multi_field.py
@@ -67,7 +67,7 @@ class MultiField(object):
         dtype = MultiField.build_dtype(dtype, domain)
         return MultiField({key: Field.from_random(random_type, domain[key],
                                                   dtype[key], **kwargs)
-                           for key in domain.keys()})
+                           for key in sorted(domain.keys())})
 
     def fill(self, fill_value):
         """Fill `self` uniformly with `fill_value`
diff --git a/nifty5/operators/exp_transform.py b/nifty5/operators/exp_transform.py
index e6e635cea2c400b0f78932ed35ac6875a6a63ddb..d9c846519b9e7e821bfb15566cba1760a5d75786 100644
--- a/nifty5/operators/exp_transform.py
+++ b/nifty5/operators/exp_transform.py
@@ -67,7 +67,7 @@ class ExpTransform(LinearOperator):
 
     def apply(self, x, mode):
         self._check_input(x, mode)
-        x = x.val
+        x = x.to_global_data()
         ndim = len(self.target.shape)
         idx = ()
         for d in range(ndim):
@@ -90,7 +90,7 @@ class ExpTransform(LinearOperator):
 
             x = xnew
             idx = (slice(None),) + idx
-        return Field(self._tgt(mode), val=x)
+        return Field.from_global_data(self._tgt(mode), x)
 
     @property
     def capability(self):
diff --git a/nifty5/operators/qht_operator.py b/nifty5/operators/qht_operator.py
index e022b990b0fc6ca2929d20ae715a6df24204c643..88da799f79c134a14c24d22b57995154884291e7 100644
--- a/nifty5/operators/qht_operator.py
+++ b/nifty5/operators/qht_operator.py
@@ -1,5 +1,6 @@
 from ..domain_tuple import DomainTuple
 from ..field import Field
+from .. import dobj
 from ..utilities import hartley
 from .linear_operator import LinearOperator
 
@@ -34,11 +35,13 @@ class QHTOperator(LinearOperator):
         x = x.val * self.domain[0].scalar_dvol()
         n = len(self.domain[0].shape)
         rng = range(n) if mode == self.TIMES else reversed(range(n))
+        # MR FIXME: this needs to be fixed properly for MPI
+        x = dobj.to_global_data(x)
         for i in rng:
             sl = (slice(None),)*i + (slice(1, None),)
             x[sl] = hartley(x[sl], axes=(i,))
 
-        return Field(self._tgt(mode), val=x)
+        return Field.from_global_data(self._tgt(mode), x)
 
     @property
     def capability(self):
diff --git a/nifty5/operators/slope_operator.py b/nifty5/operators/slope_operator.py
index 81174405321d8f3cf2575c25e75b487ef5ac04d4..3daf567233747cd6b95472efa7b4061b002ad96c 100644
--- a/nifty5/operators/slope_operator.py
+++ b/nifty5/operators/slope_operator.py
@@ -18,7 +18,7 @@ class SlopeOperator(LinearOperator):
         self.pos = np.zeros((self.ndim,) + self.target[0].shape)
 
         if self.ndim == 1:
-            self.pos[0] = self.target[0].get_k_length_array().val
+            self.pos[0] = self.target[0].get_k_length_array().to_global_data()
         else:
             shape = self.target[0].shape
             for i in range(self.ndim):
@@ -46,18 +46,19 @@ class SlopeOperator(LinearOperator):
 
         # Times
         if mode == self.TIMES:
-            inp = x.val
+            inp = x.to_global_data()
             res = self.sigmas[-1] * inp[-1]
             for i in range(self.ndim):
                 res += self.sigmas[i] * inp[i] * self.pos[i]
-            return Field(self.target, val=res)
+            return Field.from_global_data(self.target, res)
 
         # Adjoint times
         res = np.zeros(self.domain[0].shape)
-        res[-1] = np.sum(x.val) * self.sigmas[-1]
+        xglob = x.to_global_data()
+        res[-1] = np.sum(xglob) * self.sigmas[-1]
         for i in range(self.ndim):
-            res[i] = np.sum(self.pos[i] * x.val) * self.sigmas[i]
-        return Field(self.domain, val=res)
+            res[i] = np.sum(self.pos[i] * xglob) * self.sigmas[i]
+        return Field.from_global_data(self.domain, res)
 
     @property
     def capability(self):