Refs #39. added a special exception in singlephonon_sqe2dos and catch…

… it with a test case

multiphonon/backward/singlephonon_sqe2dos.py

@@ -4,6 +4,9 @@
 
 import numpy as np, histogram as H, histogram.hdf as hh, os
 
+class EnergyAxisMissingBinCenterAtZero(Exception): pass
+
+
 def sqe2dos(sqe, T, Ecutoff, elastic_E_cutoff, M, initdos=None):
     """ 
     Given a one-phonon sqe, compute dos
@@ -33,8 +36,10 @@ def sqe2dos(sqe, T, Ecutoff, elastic_E_cutoff, M, initdos=None):
     # create initial guess of dos
     Efull = sqe.E
     dE = sqe.E[1] - sqe.E[0]
+    assert dE>0, "Energy axis must be incremental"
     Eplus = Efull[Efull>-dE/2]
-    assert Eplus[0] < dE/1e6
+    if abs(Eplus[0]) > dE/1e6:
+        raise EnergyAxisMissingBinCenterAtZero('"0" must be one of the bin centers of the energy axis')
     Eplus[0] = 0.
     if initdos is None:
         initdos = guess_init_dos(Eplus, Ecutoff)

tests/backward/sqe2dos_TestCase.py

@@ -48,6 +48,17 @@ def test1b(self):
         return
 
 
+    def test1b2(self):
+        iqehist = hh.load(os.path.join(datadir, "V-iqe.h5"))
+        from multiphonon.sqe import interp
+        newiqe = interp(iqehist, newE = np.arange(-50.5, 50, 1.))
+        from multiphonon.backward.singlephonon_sqe2dos import EnergyAxisMissingBinCenterAtZero
+        with self.assertRaises(EnergyAxisMissingBinCenterAtZero):
+            DOS = sqe2dos.singlephonon_sqe2dos(
+                newiqe, T=300, Ecutoff=65., elastic_E_cutoff=(-20., 6.7), M=50.94)
+        return
+
+
     def test1c(self):
         iqehist = hh.load(os.path.join(datadir, "graphite-Ei_130-iqe.h5"))
         initdos = hh.load(os.path.join(datadir, "graphite-Ei_300-dos.h5"))