Merge branch 'development'

c/correlation.c

@@ -83,7 +83,7 @@ double EvaluateCorrelation (double Frequency, double _Complex Velocity[], int Nu
             }
         }
     }
-    return  creal(Correl) * TimeStep/(NumberOfData/Increment)/2.0;
+    return  creal(Correl) * TimeStep/(NumberOfData/Increment);
 }
 
 // Original (simple)

c/mem.c

@@ -49,7 +49,7 @@ static PyObject* MaximumEntropyMethod (PyObject* self, PyObject *arg, PyObject *
     double *PowerSpectrum  = (double*)PyArray_DATA(PowerSpectrum_object);
 
     //Declare variables
-    double Coefficients[NumberOfCoefficients];
+    double Coefficients[NumberOfCoefficients+1];
 
     // Transform complex data to double
     double *Velocity_r = (double *)malloc(NumberOfData * sizeof(double));
@@ -61,7 +61,7 @@ static PyObject* MaximumEntropyMethod (PyObject* self, PyObject *arg, PyObject *
     # pragma omp parallel for default(shared) private(AngularFrequency)
     for (int i=0; i < NumberOfFrequencies; i++) {
         AngularFrequency = Frequency[i]*2.0*M_PI;
-        PowerSpectrum[i] = FrequencyEvaluation(AngularFrequency*TimeStep, Coefficients, NumberOfCoefficients, MeanSquareDiscrepancy) * TimeStep / 2.0;
+        PowerSpectrum[i] = FrequencyEvaluation(AngularFrequency*TimeStep, Coefficients, NumberOfCoefficients, MeanSquareDiscrepancy) * TimeStep;
     }
 
     // Free python memory

dynaphopy/__init__.py

@@ -1,4 +1,4 @@
-__version__ = '1.14.4'
+__version__ = '1.14.5'
 
 import numpy as np
 import matplotlib.pyplot as plt
@@ -221,6 +221,8 @@ def plot_renormalized_phonon_dispersion_bands(self, plot_linewidths=False):
         eigenvectors = data['eigenvectors']
         linewidths = data['linewidths']
 
+        plt.suptitle('Renormalized phonon dispersion relations')
+
         sup_lim = pho_interface.get_renormalized_force_constants(renormalized_frequencies+linewidths/2,
                                                                  eigenvectors,
                                                                  self.dynamic.structure,
@@ -232,6 +234,7 @@ def plot_renormalized_phonon_dispersion_bands(self, plot_linewidths=False):
                                                                  symmetrize=self.parameters.symmetrize)
 
         if plot_linewidths:
+            plt.suptitle('Renormalized phonon dispersion relations and linewidths')
             renormalized_bands_s = pho_interface.obtain_phonon_dispersion_bands(self.dynamic.structure,
                                                                               self.parameters.band_ranges,
                                                                               force_constants=sup_lim,
@@ -258,7 +261,6 @@ def plot_renormalized_phonon_dispersion_bands(self, plot_linewidths=False):
         plt.xlabel('Wave vector')
         plt.xlim([0, self._bands[1][-1][-1]])
         plt.axhline(y=0, color='k', ls='dashed')
-        plt.suptitle('Renormalized phonon dispersion relations')
         handles, labels = plt.gca().get_legend_handles_labels()
         plt.legend([handles[0], handles[-1]], ['Harmonic', 'Renormalized'])
         plt.show()
@@ -597,7 +599,7 @@ def plot_power_spectrum_full(self):
         plt.show()
 
         total_integral = integrate.simps(self.get_power_spectrum_full(), x=self.get_frequency_range())
-        print ("Total Area (1/2 Kinetic energy <K>): {0} eV".format(total_integral))
+        print ("Total Area (Kinetic energy <K>): {0} eV".format(total_integral))
 
     def plot_power_spectrum_wave_vector(self):
         plt.suptitle('Projection onto wave vector')
@@ -606,7 +608,7 @@ def plot_power_spectrum_wave_vector(self):
         plt.ylabel('eV * ps')
         plt.show()
         total_integral = integrate.simps(self.get_power_spectrum_wave_vector(), x=self.get_frequency_range())
-        print ("Total Area (1/2 Kinetic energy <K>): {0} eV".format(total_integral))
+        print ("Total Area (Kinetic energy <K>): {0} eV".format(total_integral))
 
     def plot_power_spectrum_phonon(self):
         for i in range(self.get_power_spectrum_phonon().shape[1]):
@@ -709,14 +711,14 @@ def write_power_spectrum_full(self, file_name):
                                     self.get_power_spectrum_full()[None].T,
                                     file_name)
         total_integral = integrate.simps(self.get_power_spectrum_full(), x=self.get_frequency_range())
-        print ("Total Area (1/2 Kinetic energy <K>): {0} eV".format(total_integral))
+        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())
-        print ("Total Area (1/2 Kinetic energy <K>): {0} eV".format(total_integral))
+        print ("Total Area (Kinetic energy <K>): {0} eV".format(total_integral))
 
     def write_power_spectrum_phonon(self, file_name):
         reading.write_curve_to_file(self.get_frequency_range(),
@@ -995,11 +997,11 @@ def get_anisotropic_displacement_parameters(self, coordinate_type='uvrs', print_
         atomic_types_unique = [elements[i] for i in atom_type_index_unique]
 
         if print_on_screen:
-            print('Anisotropic displacement parameters ({0})'.format(coordinate_type))
+            print('Anisotropic displacement parameters ({0}) [relative to average atomic positions]'.format(coordinate_type))
             print('          U11          U22          U33          U23          U13          U12')
 
         anisotropic_displacements = []
-        for i, u_cart in enumerate(self.dynamic.get_mean_displacement_matrix()):
+        for i, u_cart in enumerate(self.dynamic.get_mean_displacement_matrix(use_average_positions=True)):
 
             cell = self.dynamic.structure.get_cell()
             cell_inv = np.linalg.inv(cell)
@@ -1029,6 +1031,14 @@ def get_anisotropic_displacement_parameters(self, coordinate_type='uvrs', print_
 
         return anisotropic_displacements
 
+    def get_average_atomic_positions(self):
+        print 'Average atomic positions'
+        positions_average = self.dynamic.average_positions(to_unit_cell=True)
+        elements = self.dynamic.structure.get_atomic_types()
+        for i, coordinate in enumerate(positions_average):
+            print '{0:2} '.format(elements[i]) + '{0:15.8f} {1:15.8f} {2:15.8f}'.format(*coordinate.real)
+
+
 # Support functions
 def vector_in_list(vector_test_list, vector_full_list):
     for vector_test in vector_test_list:

dynaphopy/analysis/coordinates.py

@@ -2,6 +2,7 @@
 import sys
 from dynaphopy.displacements import atomic_displacement
 
+
 def progress_bar(progress):
     bar_length = 30
     status = ""
@@ -24,14 +25,14 @@ def progress_bar(progress):
 
 #print(disp.relative_trajectory(cell, traj ,pos))
 
-#Not used (deprecated)
+#Not used (only for test)
 def relativize_trajectory(dynamic, memmap=False):
 
     cell = dynamic.get_supercell()
     number_of_atoms = dynamic.trajectory.shape[1]
     supercell = dynamic.get_supercell_matrix()
     position = dynamic.structure.get_positions(supercell=supercell)
-#    normalized_trajectory = np.zeros_like(dynamic.trajectory.real)
+    # normalized_trajectory = np.zeros_like(dynamic.trajectory.real)
     normalized_trajectory = dynamic.trajectory.copy()
 
     trajectory = dynamic.trajectory
@@ -41,15 +42,15 @@ def relativize_trajectory(dynamic, memmap=False):
     else:
         normalized_trajectory = dynamic.trajectory.copy()
 
-#    progress_bar(0)
+    # progress_bar(0)
 
     for i in range(number_of_atoms):
         normalized_trajectory[:, i, :] = atomic_displacement(trajectory[:, i, :], position[i], cell)
 
-   #     progress_bar(float(i+1)/number_of_atoms)
+   # progress_bar(float(i+1)/number_of_atoms)
     return normalized_trajectory
-#Not used anymore (deprecated)
 
+# Not used (only for test)
 def relativize_trajectory_py(dynamic):
 
     print('Using python rutine for calculating atomic displacements')
@@ -66,7 +67,7 @@ def relativize_trajectory_py(dynamic):
 
             difference = normalized_trajectory[i, j, :] - position[j]
 
-         #   difference_matrix = np.array(np.dot(np.linalg.inv(cell),(IniSep)),dtype=int)
+            # difference_matrix = np.array(np.dot(np.linalg.inv(cell),(IniSep)),dtype=int)
             difference_matrix = np.around(np.dot(np.linalg.inv(cell), difference), decimals=0)
             normalized_trajectory[i, j, :] -= np.dot(difference_matrix, cell.T) + position[j]
 
@@ -81,9 +82,6 @@ def trajectory_projection(dynamic, direction):
 
     supercell = dynamic.get_supercell_matrix()
     trajectory = dynamic.get_relative_trajectory()
-
-  #  print(trajectory)
-
     atom_type_index = dynamic.structure.get_atom_type_index(supercell=supercell)
     number_of_atom_types = dynamic.structure.get_number_of_atom_types()
 
@@ -93,7 +91,7 @@ def trajectory_projection(dynamic, direction):
         projection = np.array([])
         for i in range(0, trajectory.shape[1]):
             if atom_type_index[i] == j:
-     #           print('atom:', i, 'type:', atom_type_index[i])
+                # print('atom:', i, 'type:', atom_type_index[i])
                 projection = np.append(projection, np.dot(trajectory[:, i, :].real,
                                                           direction/np.linalg.norm(direction)))
         projections.append(projection)

dynaphopy/analysis/energy.py

@@ -2,7 +2,7 @@
 import matplotlib.pyplot as plt
 import scipy.stats as stats
 
-kb_boltzmann = 0.831446 # u * A^2 / ( ps^2 * K )
+kb_boltzmann = 0.831446  # u * A^2 / ( ps^2 * K )
 
 
 def boltzmann_distribution(trajectory, parameters):

dynaphopy/analysis/fitting/__init__.py

@@ -3,10 +3,8 @@
 from scipy.integrate import simps
 from dynaphopy.analysis.fitting import fitting_functions
 
-#from scipy.optimize import curve_fit, minimize_scalar, root
-
 h_planck = 4.135667662e-3  # eV/ps
-h_planck_bar = 6.58211951e-4  # eV/ps
+h_planck_bar = h_planck/(2.0*np.pi)  # eV/ps
 kb_boltzmann = 8.6173324e-5  # eV/K
 
 
@@ -34,7 +32,7 @@ def average_phonon(index, data, degeneracy):
             return np.average([data[:, j] for j in degeneracy[i]], axis=0)
 
 
-def phonon_fitting_analysis(original, test_frequencies_range, harmonic_frequencies=None,
+def phonon_fitting_analysis(original, ps_frequencies, harmonic_frequencies=None,
                             fitting_function_type=0,
                             show_plots=True,
                             use_degeneracy=True,
@@ -50,15 +48,14 @@ def phonon_fitting_analysis(original, test_frequencies_range, harmonic_frequenci
         if use_degeneracy:
             degeneracy = degenerate_sets(harmonic_frequencies)
             power_spectrum = average_phonon(i, original, degeneracy)
-
         else:
             power_spectrum = original[:, i]
 
         guess_height = np.max(power_spectrum)
-        guess_position = test_frequencies_range[np.argmax(power_spectrum)]
+        guess_position = ps_frequencies[np.argmax(power_spectrum)]
 
         Fitting_function_class = fitting_functions.fitting_functions[fitting_function_type]
-        fitting_function = Fitting_function_class(test_frequencies_range,
+        fitting_function = Fitting_function_class(ps_frequencies,
                                                   power_spectrum,
                                                   guess_height=guess_height,
                                                   guess_position=guess_position)
@@ -70,47 +67,49 @@ def phonon_fitting_analysis(original, test_frequencies_range, harmonic_frequenci
             widths.append(0)
             errors.append(0)
             dt_Q2_s.append(0)
+            print ('Warning: Fitting not successful in peak #{0}'.format(i+1))
             continue
 
-        frequency = fitting_parameters['peak_position']
+        position = fitting_parameters['peak_position']
         area = fitting_parameters['area']
         width = fitting_parameters['width']
         base_line = fitting_parameters['base_line']
         maximum = fitting_parameters['maximum']
         error = fitting_parameters['global_error']
 
-        total_integral = simps(power_spectrum, x=test_frequencies_range)
+        total_integral = simps(power_spectrum, x=ps_frequencies)
 
         # Calculated properties
-        dt_Q2_lor = 2 * 2 * area
-        dt_Q2_tot = 2 * 2 * total_integral
+        dt_Q2_lor = 2 * area
+        dt_Q2_tot = 2 * total_integral
 
         # Only within harmonic approximation
-        Q2_lor = dt_Q2_lor / pow(frequency * 2 * np.pi, 2)
-        Q2_tot = dt_Q2_tot / pow(frequency * 2 * np.pi,2)
+        Q2_lor = dt_Q2_lor / pow(position * 2 * np.pi, 2)
+        Q2_tot = dt_Q2_tot / pow(position * 2 * np.pi,2)
 
-        occupancy_lor = dt_Q2_lor / (frequency * h_planck_bar) - 0.5
-        occupancy_tot = dt_Q2_tot / (frequency * h_planck_bar) - 0.5
+        occupancy_lor = dt_Q2_lor / (position * h_planck_bar) - 0.5
+        occupancy_tot = dt_Q2_tot / (position * h_planck_bar) - 0.5
 
-    #    fit_temperature = dt_Q2_lor / kb_boltzmann # High temperature limit
-        fit_temperature = h_planck_bar * frequency / (kb_boltzmann * np.log((1.0 / occupancy_lor + 1.0)))
+        # fit_temperature = dt_Q2_lor / kb_boltzmann  # High temperature limit
+        fit_temperature = h_planck_bar * position / (kb_boltzmann * np.log((1.0 / occupancy_lor + 1.0)))
+        fit_temperature_tot = h_planck_bar * position / (kb_boltzmann * np.log((1.0 / occupancy_tot + 1.0)))
 
         #Print section
         print ('\nPeak # {0}'.format(i+1))
         print ('----------------------------------------------')
         print ('Width                      {0:15.6f} THz'.format(width))
-        print ('Position                   {0:15.6f} THz'.format(frequency))
-        print ('Area (1/2<K>) ({0:.10s}) {1:15.6f} eV'.format(fitting_function.curve_name, area))      # 1/2 Kinetic energy
-        print ('Area (1/2<K>) (Total)      {0:15.6f} eV'.format(total_integral))   # 1/2 Kinetic energy
+        print ('Position                   {0:15.6f} THz'.format(position))
+        print ('Area (<K>)    ({0:.10s}) {1:15.6f} eV'.format(fitting_function.curve_name, area))  # 1/2 Kinetic energy
+        print ('Area (<K>)    (Total)      {0:15.6f} eV'.format(total_integral))   # 1/2 Kinetic energy
         print ('<|dQ/dt|^2>                {0:15.6f} eV'.format(dt_Q2_lor))        # Kinetic energy
- #       print '<|dQ/dt|^2> (tot):        ', dt_Q2_tot, 'eV'        # Kinetic energy
- #       print '<|Q|^2> (lor):          ', Q2_lor, 'u * Angstrom^2'
- #       print '<|Q|^2> (tot):          ', Q2_tot, 'u * Angstrom^2'
+        # print '<|dQ/dt|^2> (tot):        ', dt_Q2_tot, 'eV'        # Kinetic energy
+        # print '<|Q|^2> (lor):          ', Q2_lor, 'eV' #  potential energy
+        # print '<|Q|^2> (tot):          ', Q2_tot, 'eV' #  potential energy
         if show_occupancy:
             print ('Occupation number          {0:15.6f}'.format(occupancy_lor))
-     #       print 'Occupation number(tot): ', occupancy_tot
-            print ('Fit temperature            {0:15.6f} K'.format(fit_temperature))
-     #       print 'Fit temperature (tot)   ', dt_Q2_tot / kb_bolzman, 'K'
+            # print ('Fit temperature            {0:15.6f} K'.format(fit_temperature))
+            print ('Fit temperature (Total)    {0:15.6f} K'.format(fit_temperature_tot))
+
         print ('Base line                  {0:15.6f} eV * ps'.format(base_line))
         print ('Maximum height             {0:15.6f} eV * ps'.format(maximum))
         print ('Fitting global error       {0:15.6f}'.format(error))
@@ -120,9 +119,9 @@ def phonon_fitting_analysis(original, test_frequencies_range, harmonic_frequenci
             print ('Peak asymmetry             {0:15.6f}'.format(asymmetry))
 
         if harmonic_frequencies is not None:
-            print ('Frequency shift            {0:15.6f} THz'.format(frequency - harmonic_frequencies[i]))
+            print ('Frequency shift            {0:15.6f} THz'.format(position - harmonic_frequencies[i]))
 
-        positions.append(frequency)
+        positions.append(position)
         widths.append(width)
         errors.append(error/maximum)
         dt_Q2_s.append(dt_Q2_lor)
@@ -136,18 +135,16 @@ def phonon_fitting_analysis(original, test_frequencies_range, harmonic_frequenci
             plt.title('Curve fitting')
 
             plt.suptitle('Phonon {0}'.format(i+1))
-            plt.text(frequency+width, guess_height/2, 'Width: ' + "{:10.4f}".format(width),
+            plt.text(position+width, guess_height/2, 'Width: ' + "{:10.4f}".format(width),
                      fontsize=12)
 
-            plt.plot(test_frequencies_range, power_spectrum,
+            plt.plot(ps_frequencies, power_spectrum,
                      label='Power spectrum')
 
-
-            plt.plot(test_frequencies_range, fitting_function.get_curve(test_frequencies_range),
+            plt.plot(ps_frequencies, fitting_function.get_curve(ps_frequencies),
                      label=fitting_function.curve_name,
                      linewidth=3)
 
-
 #            plt.plot(test_frequencies_range, dumped_harmonic(test_frequencies_range, *fit_params[:4]),
 #                     label='Lorentzian fit',
 #                     linewidth=3)
@@ -156,13 +153,12 @@ def phonon_fitting_analysis(original, test_frequencies_range, harmonic_frequenci
 #                     label=('As. Lorentzian fit' if asymmetric_peaks else 'Lorentizian fit'),
 #                     linewidth=3)
 
-            plt.axvline(x=frequency, color='k', ls='dashed')
+            plt.axvline(x=position, color='k', ls='dashed')
             plt.ylim(bottom=0)
-            plt.xlim([test_frequencies_range[0],test_frequencies_range[-1]])
-            plt.fill_between([frequency-width/2, frequency+width/2],0,plt.gca().get_ylim()[1],color='red', alpha='0.2')
+            plt.xlim([ps_frequencies[0], ps_frequencies[-1]])
+            plt.fill_between([position-width/2, position+width/2],0,plt.gca().get_ylim()[1],color='red', alpha='0.2')
             plt.legend()
 
-
     if show_plots:
         plt.show()