diff --git a/src/core/nonbonded_interactions/gb.hpp b/src/core/nonbonded_interactions/gb.hpp
index 5461dc27aa0..66dd6fcaeb5 100644
--- a/src/core/nonbonded_interactions/gb.hpp
+++ b/src/core/nonbonded_interactions/gb.hpp
@@ -30,6 +30,8 @@
 #include "nonbonded_interaction_data.hpp"
 #include "particle_data.hpp"
 #include "utils.hpp"
+#include "utils/math/int_pow.hpp"
+#include "utils/math/sqr.hpp"
 
 #ifdef GAY_BERNE
 
@@ -80,8 +82,8 @@ add_gb_pair_force(const Particle *const p1, const Particle *const p2,
   Koef1 = ia_params->GB_mu / E2;
   Koef2 = Sigma * Sigma * Sigma * 0.5;
 
-  X = 1 / (dist - Sigma + ia_params->GB_sig);
-  Xcut = 1 / (ia_params->GB_cut - Sigma + ia_params->GB_sig);
+  X = ia_params->GB_sig / (dist - Sigma + ia_params->GB_sig);
+  Xcut = ia_params->GB_sig / (ia_params->GB_cut - Sigma + ia_params->GB_sig);
 
   if (X < 1.25) { /* 1.25 corresponds to the interparticle penetration of 0.2
                     units of length.
@@ -163,47 +165,44 @@ add_gb_pair_force(const Particle *const p1, const Particle *const p2,
 inline double gb_pair_energy(const Particle *p1, const Particle *p2,
                              const IA_parameters *ia_params, const double d[3],
                              double dist) {
+  using Utils::int_pow;
+  using Utils::sqr;
+
   if (!(dist < ia_params->GB_cut))
     return 0.0;
 
-  double a, b, c, X, Xcut, Brack, BrackCut, u1x, u1y, u1z, u2x, u2y, u2z, E, E1,
-      E2, Sigma, Plus1, Minus1, Plus2, Minus2;
+  auto const e0 = ia_params->GB_eps;
+  auto const s0 = ia_params->GB_sig;
+  auto const chi1 = ia_params->GB_chi1;
+  auto const chi2 = ia_params->GB_chi2;
+  auto const mu = ia_params->GB_mu;
+  auto const nu = ia_params->GB_nu;
+  auto const r = Vector3d({d[0], d[1], d[2]}).normalize();
 
-  u1x = p1->r.calc_director()[0];
-  u1y = p1->r.calc_director()[1];
-  u1z = p1->r.calc_director()[2];
-  u2x = p2->r.calc_director()[0];
-  u2y = p2->r.calc_director()[1];
-  u2z = p2->r.calc_director()[2];
+  auto const ui = p1->r.calc_director();
+  auto const uj = p2->r.calc_director();
+  auto const uij = ui * uj;
+  auto const rui = r * ui;
+  auto const ruj = r * uj;
 
-  a = d[0] * u1x + d[1] * u1y + d[2] * u1z;
-  b = d[0] * u2x + d[1] * u2y + d[2] * u2z;
-  c = u1x * u2x + u1y * u2y + u1z * u2z;
+  auto const e1 = std::pow(1. - sqr(chi1 * uij), -0.5 * nu);
 
-  Plus1 = (a + b) / (1 + ia_params->GB_chi1 * c);
-  Plus2 = (a + b) / (1 + ia_params->GB_chi2 * c);
-  Minus1 = (a - b) / (1 - ia_params->GB_chi1 * c);
-  Minus2 = (a - b) / (1 - ia_params->GB_chi2 * c);
-  E1 = 1 / sqrt(1 - ia_params->GB_chi1 * ia_params->GB_chi1 * c * c);
-  E2 = 1 - 0.5 * (ia_params->GB_chi2 / dist / dist) *
-               (Plus2 * (a + b) + Minus2 * (a - b));
-  E = 4 * ia_params->GB_eps * pow(E1, ia_params->GB_nu) *
-      pow(E2, ia_params->GB_mu);
-  Sigma =
-      ia_params->GB_sig / sqrt(1 - 0.5 * (ia_params->GB_chi1 / dist / dist) *
-                                       (Plus1 * (a + b) + Minus1 * (a - b)));
+  auto const t1 = sqr(rui + ruj) / (1 + chi2 * uij);
+  auto const t2 = sqr(rui - ruj) / (1 - chi2 * uij);
+  auto const e2 = std::pow(1. - 0.5 * chi2 * (t1 + t2), mu);
+
+  auto const e = e0 * e1 * e2;
 
-  X = 1 / (dist - Sigma + ia_params->GB_sig);
-  Xcut = 1 / (ia_params->GB_cut - Sigma + ia_params->GB_sig);
+  auto const s = s0 / std::sqrt(1. - 0.5 * chi1 *
+                                         (sqr(rui + ruj) / (1 + chi1 * uij) +
+                                          sqr(rui - ruj) / (1 - chi1 * uij)));
 
-  Brack = X * X * X;
-  BrackCut = Xcut * Xcut * Xcut;
-  Brack = Brack * Brack;
-  BrackCut = BrackCut * BrackCut;
-  Brack = Brack * (Brack - 1);
-  BrackCut = BrackCut * (BrackCut - 1);
+  auto r_eff = [=](double r) { return (r - s + s0) / s0; };
+  auto E = [=](double r) {
+    return 4. * e * (int_pow<12>(1. / r) - int_pow<6>(1. / r));
+  };
 
-  return E * (Brack - BrackCut);
+  return E(r_eff(dist)) - E(r_eff(ia_params->GB_cut));
 }
 
 #endif
diff --git a/src/features.def b/src/features.def
index a2fb9fe91a1..9e6e3f7e180 100644
--- a/src/features.def
+++ b/src/features.def
@@ -92,7 +92,7 @@ LJCOS
 LJCOS2
 LJGEN_SOFTCORE                  implies LENNARD_JONES_GENERIC
 COS2
-GAY_BERNE
+GAY_BERNE                       depends EXPERIMENTAL_FEATURES
 SMOOTH_STEP
 HERTZIAN
 GAUSSIAN
diff --git a/testsuite/python/interactions_non-bonded.py b/testsuite/python/interactions_non-bonded.py
index 37a0e0fe8b8..12a2f8865e8 100644
--- a/testsuite/python/interactions_non-bonded.py
+++ b/testsuite/python/interactions_non-bonded.py
@@ -36,7 +36,7 @@ class InteractionsNonBondedTest(ut.TestCase):
     def setUp(self):
 
         self.system.box_l = [self.box_l] * 3
-        self.system.cell_system.skin = 0.4
+        self.system.cell_system.skin = 0.
         self.system.time_step = .1
 
         self.system.part.add(id=0, pos=self.start_pos, type=0)
@@ -585,6 +585,36 @@ def test_gaussian(self):
 
         self.system.non_bonded_inter[0, 0].gaussian.set_params(eps=0.)
 
+    @ut.skipIf(not espressomd.has_features("GAY_BERNE"), "skipped for lack of Gay Berne")
+    def test_gb(self):
+        k_1 = 1.2
+        k_2 = 2.4
+        mu = 2.
+        nu = 5.
+        sigma_0 = 1.2
+        epsilon_0 = 0.8
+        cut = 3.3 
+
+        self.system.part[:].pos = ((1, 2, 3), (2.2, 2.1, 2.9))
+        self.system.non_bonded_inter[0, 0].gay_berne.set_params(
+            sig=sigma_0, cut=cut, eps=epsilon_0, k1=k_1, k2=k_2, mu=mu, nu=nu)
+        p1 = self.system.part[0]
+        p2 = self.system.part[1]
+        p1.rotate(axis=(1, 2, 3), angle=0.3)
+        p1.rotate(axis=(1, -2, -4), angle=1.2)
+
+        r = self.system.distance_vec(p1, p2)
+
+        r_cut = r * cut / numpy.linalg.norm(r)
+        self.assertAlmostEqual(self.system.analysis.energy()["non_bonded"],
+                               tests_common.gay_berne_potential(r, p1.director, p2.director, epsilon_0, sigma_0, mu, nu, k_1, k_2) -
+                               tests_common.gay_berne_potential(r_cut, p1.director, p2.director, epsilon_0, sigma_0, mu, nu, k_1, k_2), delta=1E-14)
+
+        self.system.integrator.run(0)
+        self.system.non_bonded_inter[0, 0].gay_berne.set_params(
+            sig=sigma_0, cut=0, eps=0, k1=k_1, k2=k_2, mu=mu, nu=nu)
+        self.system.integrator.run(0)
+        self.assertEqual(self.system.analysis.energy()["non_bonded"], 0.0)
 
 if __name__ == '__main__':
     print("Features: ", espressomd.features())
diff --git a/testsuite/python/tests_common.py b/testsuite/python/tests_common.py
index c6a3058756a..0cf3e47cfca 100644
--- a/testsuite/python/tests_common.py
+++ b/testsuite/python/tests_common.py
@@ -583,6 +583,32 @@ def gaussian_force(r, eps, sig, cutoff):
     return f
 
 
+def gay_berne_potential(r_ij, u_i, u_j, epsilon_0, sigma_0, mu, nu, k_1, k_2):
+    r_normed = r_ij / np.linalg.norm(r_ij)
+    r_u_i = np.dot(r_normed, u_i)
+    r_u_j = np.dot(r_normed, u_j)
+    u_i_u_j = np.dot(u_i, u_j)
+
+    chi = (k_1**2 - 1.) / (k_1**2 + 1.)
+    chi_d = (k_2**(1. / mu) - 1) / (k_2**(1. / mu) + 1)
+
+    sigma = sigma_0 \
+        / np.sqrt(
+            (1 - 0.5 * chi * (
+             (r_u_i + r_u_j)**2 / (1 + chi * u_i_u_j) +
+             (r_u_i - r_u_j)**2 / (1 - chi * u_i_u_j))))
+
+    epsilon = epsilon_0 *\
+        (1 - chi**2 * u_i_u_j**2)**(-nu / 2.) *\
+        (1 - chi_d / 2. * (
+         (r_u_i + r_u_j)**2 / (1 + chi_d * u_i_u_j) +
+         (r_u_i - r_u_j)**2 / (1 - chi_d * u_i_u_j)))**mu
+
+    rr = np.linalg.norm((np.linalg.norm(r_ij) - sigma + sigma_0) / sigma_0)
+
+    return 4. * epsilon * (rr**-12 - rr**-6)
+
+
 class DynamicDict(dict):
 
     def __getitem__(self, key):