add compatibility with simpson scipy function

dynaphopy/__init__.py

@@ -13,7 +13,10 @@
 
 
 from dynaphopy.power_spectrum import power_spectrum_functions
-from scipy import integrate
+try:
+    from scipy.integrate import simps
+except ImportError:
+    from scipy.integrate import simpson as simps
 
 
 class Quasiparticle:
@@ -912,7 +915,7 @@ def plot_power_spectrum_full(self):
         # plt.legend()
         plt.show()
 
-        total_integral = integrate.simps(self.get_power_spectrum_full(), x=self.get_frequency_range())
+        total_integral = simps(self.get_power_spectrum_full(), x=self.get_frequency_range())
         print("Total Area (Kinetic energy <K>): {0} eV".format(total_integral))
 
     def plot_power_spectrum_wave_vector(self):
@@ -922,7 +925,7 @@ def plot_power_spectrum_wave_vector(self):
         plt.ylabel('eV * ps')
         plt.axhline(y=0, color='k', ls='dashed')
         plt.show()
-        total_integral = integrate.simps(self.get_power_spectrum_wave_vector(), x=self.get_frequency_range())
+        total_integral = simps(self.get_power_spectrum_wave_vector(), x=self.get_frequency_range())
         print("Total Area (Kinetic energy <K>): {0} eV".format(total_integral))
 
     def plot_power_spectrum_phonon(self):
@@ -1033,14 +1036,14 @@ def write_power_spectrum_full(self, file_name):
         reading.write_curve_to_file(self.get_frequency_range(),
                                     self.get_power_spectrum_full()[None].T,
                                     file_name)
-        total_integral = integrate.simps(self.get_power_spectrum_full(), x=self.get_frequency_range())
+        total_integral = simps(self.get_power_spectrum_full(), x=self.get_frequency_range())
         print("Total Area (Kinetic energy <K>): {0} eV".format(total_integral))
 
     def write_power_spectrum_wave_vector(self, file_name):
         reading.write_curve_to_file(self.get_frequency_range(),
                                     self.get_power_spectrum_wave_vector()[None].T,
                                     file_name)
-        total_integral = integrate.simps(self.get_power_spectrum_wave_vector(), x=self.get_frequency_range())
+        total_integral = simps(self.get_power_spectrum_wave_vector(), x=self.get_frequency_range())
         print("Total Area (Kinetic energy <K>): {0} eV".format(total_integral))
 
     def write_power_spectrum_phonon(self, file_name):
@@ -1260,7 +1263,7 @@ def get_thermal_properties(self, force_constants=None, normalize_dos=False):
                                                        projected_on_atom=self.parameters.project_on_atom,
                                                        NAC=self.parameters.use_NAC)
 
-        integration = integrate.simps(phonopy_dos[1], x=phonopy_dos[0]) / (
+        integration = simps(phonopy_dos[1], x=phonopy_dos[0]) / (
                 self.dynamic.structure.get_number_of_atoms() *
                 self.dynamic.structure.get_number_of_dimensions())
 
@@ -1333,7 +1336,7 @@ def display_thermal_properties(self, from_power_spectrum=False, normalize_dos=Fa
 
             power_spectrum_dos = thm.get_dos(temperature, frequency_range, self.get_power_spectrum_full(),
                                              normalization)
-            integration = integrate.simps(power_spectrum_dos, x=frequency_range) / (
+            integration = simps(power_spectrum_dos, x=frequency_range) / (
                 self.dynamic.structure.get_number_of_atoms() *
                 self.dynamic.structure.get_number_of_dimensions())
 

dynaphopy/analysis/thermal_properties.py

@@ -1,7 +1,12 @@
 import numpy as np
 import matplotlib.pylab as pl
 import warnings
-from scipy import integrate
+try:
+    from scipy.integrate import simps
+except ImportError:
+    from scipy.integrate import simpson as simps
+
+
 
 N_a = 6.022140857e23
 k_b = 1.38064852e-23  # J / K
@@ -37,7 +42,7 @@ def n(temp, freq):
     total_energy = np.nan_to_num([dos[i] * h_bar * freq * (0.5 + n(temperature, freq))
                                  for i, freq in enumerate(frequency)])
 
-    total_energy = integrate.simps(total_energy, frequency) * N_a / 1000  # KJ/K/mol
+    total_energy = simps(total_energy, frequency) * N_a / 1000  # KJ/K/mol
     return total_energy
 
 
@@ -47,7 +52,7 @@ def get_free_energy(temperature, frequency, dos):
                                  for i, freq in enumerate(frequency)])
 
     free_energy[0] = 0
-    free_energy = integrate.simps(free_energy, frequency) * N_a / 1000  # KJ/K/mol
+    free_energy = simps(free_energy, frequency) * N_a / 1000  # KJ/K/mol
     return free_energy
 
 
@@ -59,7 +64,7 @@ def n(temp, freq):
     free_energy_c = np.nan_to_num([dos[i] * -h_bar/2 *shift*(n(temperature, freq) + 1 / 2.)
                                    for i, freq in enumerate(frequency)])
 
-    free_energy_c = integrate.simps(free_energy_c, frequency) * N_a / 1000 # KJ/K/mol
+    free_energy_c = simps(free_energy_c, frequency) * N_a / 1000 # KJ/K/mol
     return free_energy_c
 
 
@@ -76,7 +81,7 @@ def n(temp, freq):
 
     free_energy_c = free_energy_1 - free_energy_2
 
-    free_energy_c = integrate.simps(free_energy_c, frequency) * N_a / 1000 # KJ/K/mol
+    free_energy_c = simps(free_energy_c, frequency) * N_a / 1000 # KJ/K/mol
     return free_energy_c
 
 
@@ -88,7 +93,7 @@ def coth(x):
     entropy = np.nan_to_num([dos[i]*(1.0 / (2. * temperature) * h_bar * freq * coth(h_bar * freq / (2 * k_b * temperature))
                                      - k_b * np.log(2 * np.sinh(h_bar * freq / (2 * k_b * temperature))))
                              for i, freq in enumerate(frequency)])
-    entropy = integrate.simps(entropy, frequency) * N_a # J/K/mol
+    entropy = simps(entropy, frequency) * N_a # J/K/mol
     return entropy
 
 # Alternative way to calculate entropy (not used)
@@ -100,7 +105,7 @@ def n(temp, freq):
     entropy = np.nan_to_num([dos[i] * k_b * ((n(temperature, freq) + 1) * np.log(n(temperature, freq) + 1)
                                              - n(temperature, freq) * np.log(n(temperature, freq)))
                          for i, freq in enumerate(frequency)])
-    entropy = integrate.simps(entropy, frequency) * N_a # J/K/mol
+    entropy = simps(entropy, frequency) * N_a # J/K/mol
     return entropy
 
 
@@ -111,7 +116,7 @@ def z(temp, freq):
 
     c_v = np.nan_to_num([dos[i] * k_b * pow(z(temperature, freq), 2) * np.exp(z(temperature, freq)) / pow(np.exp(z(temperature, freq)) - 1, 2)
                          for i, freq in enumerate(frequency)])
-    c_v = integrate.simps(c_v, frequency) * N_a # J/K/mol
+    c_v = simps(c_v, frequency) * N_a # J/K/mol
 
     return c_v
 
@@ -171,7 +176,7 @@ def z(temp, freq):
     print ('Entropy: {0} J/K/mol'.format(entropy))
     print ('Cv: {0} J/K/mol'.format(c_v))
     print (np.trapz(dos_r, x=frequency_r))/(8*3)
-    print (integrate.simps(dos_r,x=frequency_r)/(8*3))
+    print (simps(dos_r,x=frequency_r)/(8*3))
 
     print ('\nFrom MD')
     print ('-------------------------')
@@ -183,7 +188,7 @@ def z(temp, freq):
     print ('Entropy: {0} J/K/mol'.format(entropy))
     print ('Cv: {0} J/K/mol'.format(c_v))
     print (np.trapz(power_spectrum, x=frequency_p))/(8*3)
-    print (integrate.simps(power_spectrum, x=frequency_p))/(8*3)
+    print (simps(power_spectrum, x=frequency_p))/(8*3)
 
     print ('\nHARMONIC')
     print ('-------------------------')
@@ -195,4 +200,4 @@ def z(temp, freq):
     print ('Entropy: {0} J/K/mol'.format(entropy))
     print ('Cv: {0} J/K/mol'.format(c_v))
     print (np.trapz(dos, x=frequency)/(8*3))
-    print (integrate.simps(dos, x=frequency)/(8*3))
+    print (simps(dos, x=frequency)/(8*3))