It seems the problem was not with the primitive estimator but when initializing all particles at 0 0 0.

When initialized randomly it works.

src/REPLICA/fix_pimdb_langevin.cpp

@@ -132,8 +132,6 @@ void FixPIMDBLangevin::compute_t_prim()
 {
     if (universe->iworld == 0)
         t_prim = bosonic_exchange.prim_estimator();
-    else
-        t_prim = -bosonic_exchange.get_total_spring_energy_for_bead();
 }
 
 /* ---------------------------------------------------------------------- */