diff --git a/bin/matdyn2fqha.py b/bin/matdyn2fqha.py
index d8d91c59..4823e9a2 100755
--- a/bin/matdyn2fqha.py
+++ b/bin/matdyn2fqha.py
@@ -34,14 +34,14 @@
 
 import sys
 import numpy as np
-from scipy.integrate import simps
+from scipy.integrate import simpson as simps
 
 filename = sys.argv[1]
 arr = np.loadtxt(filename)
 freq = arr[:,0]
 dos = arr[:,1]
 
-integral = simps(dos, freq)
+integral = simps(dos, x=freq)
 
 natom = integral / 3.
 nstep = len(freq)
diff --git a/src/pwtools/num.py b/src/pwtools/num.py
index ea04ddb1..b5963110 100644
--- a/src/pwtools/num.py
+++ b/src/pwtools/num.py
@@ -3,7 +3,7 @@
 from math import sqrt, sin, cos, radians, pi
 import scipy.optimize as optimize
 from scipy.interpolate import bisplrep, bisplev, splev, splrep
-from scipy.integrate import simps
+from scipy.integrate import simpson as simps
 from pwtools import _flib
 import warnings
 
@@ -87,7 +87,7 @@ def norm_int(y, x, area=1.0, scale=True, func=simps):
         different scales.
     func : callable
         Function to do integration (like scipy.integrate.{simps,trapz,...}
-        Called as ``func(y,x)``. Default: simps
+        Called as ``func(y,x=x)``. Default: simps
 
     Returns
     -------
@@ -107,7 +107,7 @@ def norm_int(y, x, area=1.0, scale=True, func=simps):
         fx = fy = 1.0
         sx, sy = x, y
     # Area under unscaled y(x).
-    _area = func(sy, sx) * fx * fy
+    _area = func(sy, x=sx) * fx * fy
     return y * area / _area
 
 
@@ -1097,7 +1097,7 @@ def a2_to_an(self):
             a = np.unique(self.a2[:, colidx])
             axes.append(a)
             dims.append(len(a))
-        assert np.product(dims) == self.a2.shape[0]
+        assert np.prod(dims) == self.a2.shape[0]
         idx = itertools.product(*tuple(map(range, dims)))
         an = np.empty(dims, dtype=self.a2.dtype)
         # an[1,2,3] == an[(1,2,3)], need way to eliminate loop over index array
diff --git a/src/pwtools/thermo.py b/src/pwtools/thermo.py
index f625ba97..61a544c8 100644
--- a/src/pwtools/thermo.py
+++ b/src/pwtools/thermo.py
@@ -3,7 +3,7 @@
 import warnings
 
 import numpy as np
-from scipy.integrate import simps, trapz
+from scipy.integrate import trapezoid as trapz
 from pwtools.constants import kb, hplanck, R, pi, c0, Ry_to_J, eV,\
     eV_by_Ang3_to_GPa
 from pwtools.verbose import verbose
@@ -48,7 +48,7 @@ def __init__(self, freq, dos, T=None, temp=None, skipfreq=False,
             If not None, then re-normalize the area int(freq) dos to `dosarea`,
             after `skipfreq` was applied if used.
         integrator : callable
-            Function which integrates x-y data. Called as ``integrator(y,x)``,
+            Function which integrates x-y data. Called as ``integrator(y,x=x)``,
             like ``scipy.integrate.{trapz,simps}``. Usually, `trapz` is
             numerically more stable for weird DOS data and accurate enough if
             the frequency axis resolution is good.
@@ -158,7 +158,7 @@ def _norm_int(self, y, x, area):
         fy = np.abs(y).max()
         sx = x / fx
         sy = y / fy
-        _area = self.integrator(sy, sx) * fx * fy
+        _area = self.integrator(sy, x=sx) * fx * fy
         return y*area/_area
 
     def _printwarn(self, msg):
diff --git a/test/test_norm_int.py b/test/test_norm_int.py
index b3ebd3dd..61b8c78b 100644
--- a/test/test_norm_int.py
+++ b/test/test_norm_int.py
@@ -1,5 +1,5 @@
 import numpy as np
-from scipy.integrate import simps
+from scipy.integrate import simpson as simps
 from pwtools.num import norm_int
 
 def test_norm_int():
@@ -9,4 +9,4 @@ def test_norm_int():
 
     for scale in [True, False]:
         yy = norm_int(y, x, area=10.0, scale=scale)
-        assert np.allclose(simps(yy,x), 10.0)
+        assert np.allclose(simps(yy,x=x), 10.0)
diff --git a/test/test_qha.py b/test/test_qha.py
index ed2371c5..b6227c71 100644
--- a/test/test_qha.py
+++ b/test/test_qha.py
@@ -7,7 +7,7 @@
 import numpy as np
 from numpy.testing import assert_array_almost_equal as aaae
 
-from scipy.integrate import simps, trapz
+from scipy.integrate import simpson as simps
 from pwtools.thermo import HarmonicThermo
 from pwtools import common
 from pwtools.constants import Ry_to_J, eV, Ry, kb
@@ -91,4 +91,4 @@ def test_qha():
     area = np.random.rand()*10
     ha = HarmonicThermo(pdos[:,0], pdos[:,1], skipfreq=True, dosarea=area,
                         integrator=simps)
-    assert np.allclose(simps(ha.dos, ha.f), area)
+    assert np.allclose(simps(ha.dos, x=ha.f), area)
diff --git a/test/test_trajectory.py b/test/test_trajectory.py
index 2f74992b..b9eaaeb9 100644
--- a/test/test_trajectory.py
+++ b/test/test_trajectory.py
@@ -4,7 +4,7 @@
 import copy
 
 import numpy as np
-from scipy.signal import hann
+from scipy.signal.windows import hann
 
 from pwtools.crys import Trajectory, Structure
 from pwtools import crys, constants