src/equilibrium.f90: fix floating point exception

If frhopol is called with a very small ρ_p, less than then step size of
ε_machine^⅓, the two sided finite step differentiation evaluates frhotor
at ρ_p < 0, producing a floating point exception.
In this case we use the single- sided definition.

src/equilibrium.f90

@@ -1155,7 +1155,12 @@ contains
           f(1) = frhotor(x(1)) - rho_t
         else
           ! return f'(x), computed numerically
+          if (x(1) - e > 0) then
             df(1,1) = (frhotor(x(1) + e) - frhotor(x(1) - e)) / (2*e)
+          else
+            ! single-sided when close to ρ=0
+            df(1,1) = (frhotor(x(1) + e) - frhotor(x(1))) / e
+          end if
         end if
       end subroutine