diff --git a/c/correlation.c b/c/correlation.c index 2286cef..2ee82fc 100644 --- a/c/correlation.c +++ b/c/correlation.c @@ -109,8 +109,8 @@ static PyObject* correlation_par (PyObject* self, PyObject *arg, PyObject *keywo static char *kwlist[] = {"frequency", "velocity", "time_step", "step", "integration_method", NULL}; if (!PyArg_ParseTupleAndKeywords(arg, keywords, "OOd|ii", kwlist, &frequency_obj, &velocity_obj, &TimeStep, &Increment, &IntMethod)) return NULL; - PyObject *velocity_array = PyArray_FROM_OTF(velocity_obj, NPY_CDOUBLE, NPY_IN_ARRAY); - PyObject *frequency_array = PyArray_FROM_OTF(frequency_obj, NPY_DOUBLE, NPY_IN_ARRAY); + PyObject *velocity_array = PyArray_FROM_OTF(velocity_obj, NPY_CDOUBLE, NPY_ARRAY_IN_ARRAY); + PyObject *frequency_array = PyArray_FROM_OTF(frequency_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); if (velocity_array == NULL || frequency_array == NULL ) { Py_XDECREF(velocity_array); @@ -118,17 +118,16 @@ static PyObject* correlation_par (PyObject* self, PyObject *arg, PyObject *keywo return NULL; } - _Dcomplex *Velocity = (_Dcomplex *)PyArray_DATA(velocity_array); - double *Frequency = (double*)PyArray_DATA(frequency_array); - int NumberOfData = (int)PyArray_DIM(velocity_array, 0); - int NumberOfFrequencies = (int)PyArray_DIM(frequency_array, 0); + _Dcomplex *Velocity = (_Dcomplex *)PyArray_DATA((PyArrayObject *)velocity_array); + double *Frequency = (double*)PyArray_DATA((PyArrayObject *)frequency_array); + npy_intp NumberOfData = PyArray_DIM((PyArrayObject *)velocity_array, 0); + npy_intp NumberOfFrequencies = PyArray_DIM((PyArrayObject *)frequency_array, 0); - - //Create new numpy array for storing result + // Create new numpy array for storing result PyArrayObject *PowerSpectrum_object; - int dims[1]={NumberOfFrequencies}; - PowerSpectrum_object = (PyArrayObject *) PyArray_FromDims(1,dims,NPY_DOUBLE); - double *PowerSpectrum = (double*)PyArray_DATA(PowerSpectrum_object); + npy_intp dims[] = {NumberOfFrequencies}; + PowerSpectrum_object = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE); + double *PowerSpectrum = (double *)PyArray_DATA(PowerSpectrum_object); // Maximum Entropy Method Algorithm if (IntMethod < 0 || IntMethod > 1) { diff --git a/c/displacements.c b/c/displacements.c index 44b886b..4b19f79 100644 --- a/c/displacements.c +++ b/c/displacements.c @@ -99,9 +99,9 @@ static PyObject *atomic_displacements(PyObject *self, PyObject *arg, PyObject *k static char *kwlist[] = {"trajectory", "positions", "cell", NULL}; if (!PyArg_ParseTupleAndKeywords(arg, keywords, "OOO", kwlist, &Trajectory_obj, &Positions_obj, &Cell_obj)) return NULL; - PyObject *Trajectory_array = PyArray_FROM_OTF(Trajectory_obj, NPY_CDOUBLE, NPY_IN_ARRAY); - PyObject *Positions_array = PyArray_FROM_OTF(Positions_obj, NPY_DOUBLE, NPY_IN_ARRAY); - PyObject *Cell_array = PyArray_FROM_OTF(Cell_obj, NPY_DOUBLE, NPY_IN_ARRAY); + PyObject *Trajectory_array = PyArray_FROM_OTF(Trajectory_obj, NPY_CDOUBLE, NPY_ARRAY_IN_ARRAY); + PyObject *Positions_array = PyArray_FROM_OTF(Positions_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); + PyObject *Cell_array = PyArray_FROM_OTF(Cell_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY); if (Cell_array == NULL || Trajectory_array == NULL) { Py_XDECREF(Cell_array); @@ -112,22 +112,21 @@ static PyObject *atomic_displacements(PyObject *self, PyObject *arg, PyObject *k // double _Complex *Cell = (double _Complex*)PyArray_DATA(Cell_array); - _Dcomplex *Trajectory = (_Dcomplex*)PyArray_DATA(Trajectory_array); - double *Positions = (double*)PyArray_DATA(Positions_array); + _Dcomplex *Trajectory = (_Dcomplex*)PyArray_DATA((PyArrayObject *)Trajectory_array); + double *Positions = (double*)PyArray_DATA((PyArrayObject *)Positions_array); - int NumberOfData = (int)PyArray_DIM(Trajectory_array, 0); - int NumberOfDimensions = (int)PyArray_DIM(Cell_array, 0); + npy_intp NumberOfData = PyArray_DIM((PyArrayObject *)Trajectory_array, 0); + npy_intp NumberOfDimensions = PyArray_DIM((PyArrayObject *)Cell_array, 0); // Create new Numpy array to store the result _Dcomplex **Displacement; PyArrayObject *Displacement_object; - npy_intp dims[]={NumberOfData, NumberOfDimensions}; - Displacement_object=(PyArrayObject *) PyArray_SimpleNew(2, dims, NPY_CDOUBLE); - Displacement=pymatrix_to_c_array_complex( Displacement_object); + npy_intp dims[] = {NumberOfData, NumberOfDimensions}; + Displacement_object = (PyArrayObject *)PyArray_SimpleNew(2, dims, NPY_CDOUBLE); + Displacement = pymatrix_to_c_array_complex(Displacement_object); - -// Create a pointer array from cell matrix (to be improved) + // Create a pointer array from cell matrix (to be improved) double **Cell_c = pymatrix_to_c_array_real((PyArrayObject *) Cell_array); /* @@ -196,11 +195,12 @@ static PyObject *atomic_displacements(PyObject *self, PyObject *arg, PyObject *k static _Dcomplex **pymatrix_to_c_array_complex(PyArrayObject *array) { - long n=(*array).dimensions[0]; - long m=(*array).dimensions[1]; + npy_intp n = PyArray_DIM(array, 0); + npy_intp m = PyArray_DIM(array, 1); + _Dcomplex ** c = malloc(n*sizeof(_Dcomplex)); - _Dcomplex *a = (_Dcomplex *) (*array).data; /* pointer to array data as double _Complex */ + _Dcomplex *a = (_Dcomplex *)PyArray_DATA(array); /* pointer to array data as double _Complex */ for ( int i=0; i): {0} eV".format(total_integral)) def plot_power_spectrum_wave_vector(self): @@ -920,7 +920,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 = integrate.simpson(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): @@ -1031,14 +1031,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 = integrate.simpson(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 = integrate.simpson(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 +1260,7 @@ def get_thermal_properties(self, force_constants=None): free_energy = thm.get_free_energy(temperature, phonopy_dos[0], phonopy_dos[1]) entropy = thm.get_entropy(temperature, phonopy_dos[0], phonopy_dos[1]) c_v = thm.get_cv(temperature, phonopy_dos[0], phonopy_dos[1]) - integration = integrate.simps(phonopy_dos[1], x=phonopy_dos[0]) / ( + integration = integrate.simpson(phonopy_dos[1], x=phonopy_dos[0]) / ( self.dynamic.structure.get_number_of_atoms() * self.dynamic.structure.get_number_of_dimensions()) total_energy = thm.get_total_energy(temperature, phonopy_dos[0], phonopy_dos[1]) @@ -1319,7 +1319,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 = integrate.simpson(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/fitting/__init__.py b/dynaphopy/analysis/fitting/__init__.py index 61b1249..bcaf0ad 100644 --- a/dynaphopy/analysis/fitting/__init__.py +++ b/dynaphopy/analysis/fitting/__init__.py @@ -1,6 +1,6 @@ import numpy as np import matplotlib.pyplot as plt -from scipy.integrate import simps +from scipy.integrate import simpson from dynaphopy.analysis.fitting import fitting_functions h_planck = 4.135667662e-3 # eV/ps @@ -79,7 +79,7 @@ def phonon_fitting_analysis(original, ps_frequencies, maximum = fitting_parameters['maximum'] error = fitting_parameters['global_error'] - total_integral = simps(power_spectrum, x=ps_frequencies) + total_integral = simpson(power_spectrum, x=ps_frequencies) # Calculated properties dt_Q2_lor = 2 * area diff --git a/dynaphopy/analysis/thermal_properties.py b/dynaphopy/analysis/thermal_properties.py index 2e67079..bc62925 100644 --- a/dynaphopy/analysis/thermal_properties.py +++ b/dynaphopy/analysis/thermal_properties.py @@ -37,7 +37,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 = integrate.simpson(total_energy, x=frequency) * N_a / 1000 # KJ/K/mol return total_energy @@ -47,7 +47,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 = integrate.simpson(free_energy, x=frequency) * N_a / 1000 # KJ/K/mol return free_energy @@ -59,7 +59,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 = integrate.simpson(free_energy_c, x=frequency) * N_a / 1000 # KJ/K/mol return free_energy_c @@ -76,7 +76,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 = integrate.simpson(free_energy_c, x=frequency) * N_a / 1000 # KJ/K/mol return free_energy_c @@ -88,7 +88,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 = integrate.simpson(entropy, x=frequency) * N_a # J/K/mol return entropy # Alternative way to calculate entropy (not used) @@ -100,7 +100,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 = integrate.simpson(entropy, x=frequency) * N_a # J/K/mol return entropy @@ -111,7 +111,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 = integrate.simpson(c_v, x=frequency) * N_a # J/K/mol return c_v @@ -171,7 +171,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 (integrate.simpson(dos_r,x=frequency_r)/(8*3)) print ('\nFrom MD') print ('-------------------------') @@ -183,7 +183,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 (integrate.simpson(power_spectrum, x=frequency_p))/(8*3) print ('\nHARMONIC') print ('-------------------------') @@ -195,4 +195,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 (integrate.simpson(dos, x=frequency)/(8*3)) \ No newline at end of file diff --git a/dynaphopy/atoms.py b/dynaphopy/atoms.py index de746e5..3d5b51e 100644 --- a/dynaphopy/atoms.py +++ b/dynaphopy/atoms.py @@ -187,7 +187,7 @@ def get_force_sets(self): if not isinstance(self._force_sets,type(None)): force_atoms_file = self._force_sets.get_dict()['natom'] - force_atoms_input = np.product(np.diagonal(self._force_sets.get_supercell())) * self.get_number_of_atoms() + force_atoms_input = np.prod(np.diagonal(self._force_sets.get_supercell())) * self.get_number_of_atoms() if force_atoms_file != force_atoms_input: print("Error: FORCE_SETS file does not match with SUPERCELL MATRIX") diff --git a/dynaphopy/dynamics.py b/dynaphopy/dynamics.py index 61b5103..9771bcf 100644 --- a/dynaphopy/dynamics.py +++ b/dynaphopy/dynamics.py @@ -108,7 +108,7 @@ def crop_trajectory(self, last_steps): def get_number_of_atoms(self): if self._number_of_atoms is None: - self._number_of_atoms = self.structure.get_number_of_atoms()*np.product(self.get_supercell_matrix()) + self._number_of_atoms = self.structure.get_number_of_atoms()*np.prod(self.get_supercell_matrix()) return self._number_of_atoms def set_time(self, time): diff --git a/setup.py b/setup.py index 9296cbb..4dae0c8 100644 --- a/setup.py +++ b/setup.py @@ -77,7 +77,7 @@ def get_version_number(): 'scripts/fitdata', 'scripts/qha_extract', 'scripts/rfc_calc'], - install_requires=['phonopy', 'numpy', 'scipy', 'matplotlib'] + ["windows-curses"] if sys.platform in ["win32", "cygwin"] else [], + install_requires=['phonopy', 'numpy', 'scipy', 'matplotlib'] + (["windows-curses"] if sys.platform in ["win32", "cygwin"] else []), license='MIT License', long_description=open('README.md').read(), long_description_content_type='text/markdown',