diff --git a/nifty6/operators/energy_operators.py b/nifty6/operators/energy_operators.py
index c715fc90e8fe92d46df2b689da5d6faf4098691b..597d66bab2f81f3bca0d31c7a788de45d6dfcfb7 100644
--- a/nifty6/operators/energy_operators.py
+++ b/nifty6/operators/energy_operators.py
@@ -138,7 +138,7 @@ class VariableCovarianceGaussianEnergy(EnergyOperator):
             return res
         mf = {self._r: x.val[self._icov], self._icov: .5*x.val[self._icov]**(-2)}
         metric = makeOp(MultiField.from_dict(mf))
-        return res.add_metric(SandwichOperator(x.jac, metric))
+        return res.add_metric(SandwichOperator.make(x.jac, metric))
 
 
 class GaussianEnergy(EnergyOperator):
diff --git a/test/test_energy_gradients.py b/test/test_energy_gradients.py
index adea87054c2187ddbb325733bceb15f0afe6b019..1b0972c22b7a618e5e255ceeb68087046eee7972 100644
--- a/test/test_energy_gradients.py
+++ b/test/test_energy_gradients.py
@@ -48,6 +48,7 @@ def test_variablecovariancegaussian(field):
     mf = ift.MultiField.from_dict(dc)
     energy = ift.VariableCovarianceGaussianEnergy(field.domain, 'a', 'b')
     ift.extra.check_jacobian_consistency(energy, mf, tol=1e-6)
+    energy(ift.Linearization.make_var(mf, want_metric=True)).metric.draw_sample()
 
 
 def test_gaussian(field):