replace stability error with warning message in Lorentzian susceptibi…

…lities (NanoComp#150)

* fix bug due to rounding error in lorentzian_unstable

* replace the error with a warning message instead

* rm test

src/susceptibility.cpp

@@ -152,7 +152,9 @@ void *lorentzian_susceptibility::copy_internal_data(void *data) const {
    Note that the pole satisfies the quadratic equation:
             (z + 1/z - 2)/dt^2 + g*(z - 1/z)/(2*dt) + w^2 = 0
    where w = 2*pi*omega_0 and g = 2*pi*gamma.   It is just a little
-   algebra from this to get the condition for a root with |z| > 1. */
+   algebra from this to get the condition for a root with |z| > 1.
+   
+   FIXME: this test seems to be too conservative (issue #12) */
 static bool lorentzian_unstable(double omega_0, double gamma, double dt) {
   double w = 2*pi*omega_0, g = 2*pi*gamma;
   double g2 = g*dt/2, w2 = (w*dt)*(w*dt);
@@ -177,10 +179,8 @@ void lorentzian_susceptibility::update_P
   const double omega0dtsqr_denom = no_omega_0_denominator ? 0 : omega0dtsqr;
   (void) W_prev; // unused;
 
-  if (!no_omega_0_denominator && gamma >= 0
-      && lorentzian_unstable(omega_0, gamma, dt))
-    abort("Lorentzian pole at too high a frequency %g for stability with dt = %g: reduce the Courant factor, increase the resolution, or use a different dielectric model\n", omega_0, dt);
-
+  // TODO: add back lorentzian_unstable(omega_0, gamma, dt) if we can improve the stability test
+
   FOR_COMPONENTS(c) DOCMP2 if (d->P[c][cmp]) {
     const realnum *w = W[c][cmp], *s = sigma[c][component_direction(c)];
     if (w && s) {