From 43ec15080295bdaed5a3d151f4da60e1f84ad1d7 Mon Sep 17 00:00:00 2001 From: Abel Carreras Date: Wed, 6 Nov 2024 15:22:24 +0100 Subject: [PATCH] add compatibility with simpson scipy function --- dynaphopy/__init__.py | 17 +++++++++------ dynaphopy/analysis/thermal_properties.py | 27 ++++++++++++++---------- 2 files changed, 26 insertions(+), 18 deletions(-) diff --git a/dynaphopy/__init__.py b/dynaphopy/__init__.py index 4792b1a..b3e9189 100644 --- a/dynaphopy/__init__.py +++ b/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 ): {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 ): {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 ): {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 ): {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()) diff --git a/dynaphopy/analysis/thermal_properties.py b/dynaphopy/analysis/thermal_properties.py index 2e67079..4e1eaa2 100644 --- a/dynaphopy/analysis/thermal_properties.py +++ b/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)) \ No newline at end of file + print (simps(dos, x=frequency)/(8*3)) \ No newline at end of file