diff --git a/nifty6/library/correlated_fields.py b/nifty6/library/correlated_fields.py
index 4f5932240159997f8bdd051757ac36585f1c4e38..9955bb1b7253dd6ca3b7ea3fbab9ad297886482b 100644
--- a/nifty6/library/correlated_fields.py
+++ b/nifty6/library/correlated_fields.py
@@ -18,6 +18,7 @@
 
 import numpy as np
 
+from ..logger import logger
 from ..domain_tuple import DomainTuple
 from ..domains.power_space import PowerSpace
 from ..domains.unstructured_domain import UnstructuredDomain
@@ -501,7 +502,7 @@ class CorrelatedFieldMaker:
             mean = sc.mean.to_global_data()
             stddev = sc.var.sqrt().to_global_data()
             for m, s in zip(mean.flatten(), stddev.flatten()):
-                print('{}: {:.02E} ± {:.02E}'.format(kk, m, s))
+                logger.info('{}: {:.02E} ± {:.02E}'.format(kk, m, s))
 
     def moment_slice_to_average(self, fluctuations_slice_mean, nsamples=1000):
         fluctuations_slice_mean = float(fluctuations_slice_mean)
diff --git a/nifty6/sugar.py b/nifty6/sugar.py
index c60970bca17ab412d2e6a3be4062cf8009c744bb..21bbfa28ecfaca5fc26afa5417ae23b926d24fe3 100644
--- a/nifty6/sugar.py
+++ b/nifty6/sugar.py
@@ -20,6 +20,7 @@ from time import time
 
 import numpy as np
 
+from .logger import logger
 from . import dobj, utilities
 from .domain_tuple import DomainTuple
 from .domains.power_space import PowerSpace
@@ -462,46 +463,46 @@ def exec_time(obj, want_metric=True):
     if isinstance(obj, Energy):
         t0 = time()
         obj.at(0.99*obj.position)
-        print('Energy.at():', time() - t0)
+        logger.info('Energy.at(): {}'.format(time() - t0))
 
         t0 = time()
         obj.value
-        print('Energy.value:', time() - t0)
+        logger.info('Energy.value: {}'.format(time() - t0))
         t0 = time()
         obj.gradient
-        print('Energy.gradient:', time() - t0)
+        logger.info('Energy.gradient: {}'.format(time() - t0))
         t0 = time()
         obj.metric
-        print('Energy.metric:', time() - t0)
+        logger.info('Energy.metric: {}'.format(time() - t0))
 
         t0 = time()
         obj.apply_metric(obj.position)
-        print('Energy.apply_metric:', time() - t0)
+        logger.info('Energy.apply_metric: {}'.format(time() - t0))
 
         t0 = time()
         obj.metric(obj.position)
-        print('Energy.metric(position):', time() - t0)
+        logger.info('Energy.metric(position): {}'.format(time() - t0))
     elif isinstance(obj, Operator):
         want_metric = bool(want_metric)
         pos = from_random('normal', obj.domain)
         t0 = time()
         obj(pos)
-        print('Operator call with field:', time() - t0)
+        logger.info('Operator call with field: {}'.format(time() - t0))
 
         lin = Linearization.make_var(pos, want_metric=want_metric)
         t0 = time()
         res = obj(lin)
-        print('Operator call with linearization:', time() - t0)
+        logger.info('Operator call with linearization: {}'.format(time() - t0))
 
         if isinstance(obj, EnergyOperator):
             t0 = time()
             res.gradient
-            print('Gradient evaluation:', time() - t0)
+            logger.info('Gradient evaluation: {}'.format(time() - t0))
 
             if want_metric:
                 t0 = time()
                 res.metric(pos)
-                print('Metric apply:', time() - t0)
+                logger.info('Metric apply: {}'.format(time() - t0))
     else:
         raise TypeError