Merge branch 'development'

README

@@ -1,5 +1,5 @@
 ---------------------------------------------------------
-DYNAPHOPY 1.9
+DYNAPHOPY 1.10
 ---------------------------------------------------------
 Program to calculate crystal anharmonic frequency shifts
 and phonon linewidths from molecular dynamics simulations

dynaphopy/__init__.py

@@ -1,5 +1,5 @@
+__version__='1.10'
 
-__version__='1.9'
 
 import numpy as np
 import matplotlib.pyplot as plt
@@ -22,11 +22,8 @@
     3: [power_spectrum.get_fft_fftw_spectra,           'Fast Fourier transform (FFTW)']
 }
 
-class Calculation:
 
-    @property
-    def __version__(self):
-        return '1.8'
+class Calculation:
 
     def __init__(self,
                  dynamic,
@@ -44,7 +41,6 @@ def __init__(self,
         self._bands = None
         self._renormalized_bands = None
         self._renormalized_force_constants = None
-
         self._parameters = parameters.Parameters()
         self.crop_trajectory(last_steps)
       #  print('Using {0} time steps for calculation'.format(len(self.dynamic.velocity)))
@@ -256,16 +252,18 @@ def plot_dos_phonopy(self):
         plt.suptitle('Density of states')
         plt.show()
 
-    def check_commensurate(self, q_vector, decimals=4):
+
+    def check_commensurate(self, q_point, decimals=4):
         super_cell= self.dynamic.get_super_cell_matrix()
 
         commensurate = False
         primitive_matrix = self.dynamic.structure.get_primitive_matrix()
-        q_point_unit_cell = np.dot(q_vector, np.linalg.inv(primitive_matrix))
-        q_point_unit_cell = np.multiply(q_point_unit_cell, super_cell)*2
-        q_point_unit_cell = np.around(q_point_unit_cell, decimals=decimals)
 
-        if np.all(np.equal(np.mod(q_point_unit_cell, 1), 0)):
+        transform = np.dot(q_point, np.linalg.inv(primitive_matrix))
+        transform = np.multiply(transform, super_cell)
+        transform = np.around(transform, decimals=decimals)
+
+        if np.all(np.equal(np.mod(transform, 1), 0)):
             commensurate = True
 
         return commensurate
@@ -276,7 +274,7 @@ def get_vc(self):
             print("Projecting into wave vector")
             #Check if commensurate point
             if not self.check_commensurate(self.get_reduced_q_vector()):
-                print("warning! Defined wave vector may not be a commensurate q-point in phonopy cell")
+                print("warning! Defined wave vector is not a commensurate q-point in MD supercell")
             self._vc = projection.project_onto_wave_vector(self.dynamic, self.get_q_vector())
         return self._vc
 
@@ -452,12 +450,7 @@ def plot_power_spectrum_full(self):
         handles2, labels = ax2.get_legend_handles_labels()
 
         handles = handles1 + handles2
-
-      #  print(handles)
-
-        plt.legend(handles, ['Molecular dynamics PS', 'DoS (Harmonic Aprox.)', 'DoS (renormalized)'])
-
-
+        plt.legend(handles, ['Molecular dynamics', 'DoS (Harmonic)', 'DoS (Renormalized)'])
 #        plt.legend()
         plt.show()
 
@@ -557,11 +550,15 @@ def write_power_spectrum_full(self, file_name):
         reading.write_correlation_to_file(self.get_frequency_range(),
                                           self.get_power_spectrum_direct()[None].T,
                                           file_name)
+        total_integral = np.trapz(self.get_power_spectrum_direct(), x=self.get_frequency_range())/(2 * np.pi)
+        print ("Total Area (1/2 Kinetic energy <K>): {0} eV".format(total_integral))
 
     def write_power_spectrum_wave_vector(self, file_name):
         reading.write_correlation_to_file(self.get_frequency_range(),
                                           self.get_power_spectrum_wave_vector()[None].T,
                                           file_name)
+        total_integral = np.trapz(self.get_power_spectrum_wave_vector(), x=self.get_frequency_range())/(2 * np.pi)
+        print ("Total Area (1/2 Kinetic energy <K>): {0} eV".format(total_integral))
 
     def write_power_spectrum_phonon(self,file_name):
         reading.write_correlation_to_file(self.get_frequency_range(),
@@ -617,7 +614,11 @@ def get_algorithm_list(self):
     def get_renormalized_constants(self):
 
         if self._renormalized_force_constants is None:
-            com_points = pho_interface.get_commensurate_points(self.dynamic.structure)
+            if self.parameters.use_MD_cell_commensurate:
+                self.dynamic.structure.set_super_cell_phonon_renormalized(np.diag(self.dynamic.get_super_cell_matrix()))
+
+            com_points = pho_interface.get_commensurate_points(self.dynamic.structure,
+                                                               custom_supercell=self.dynamic.structure.get_super_cell_phonon_renormalized())
 
             initial_reduced_q_point = self.get_reduced_q_vector()
 
@@ -729,50 +730,59 @@ def display_thermal_properties(self, temperature=None, from_power_spectrum=False
         c_v = thm.get_cv(temperature, phonopy_dos[0], phonopy_dos[1])
         integration = np.trapz(phonopy_dos[1], x=phonopy_dos[0])/(self.dynamic.structure.get_number_of_atoms()*
                                                        self.dynamic.structure.get_number_of_dimensions())
-        harmonic_properties = [free_energy, entropy, c_v, integration]
+        total_energy = thm.get_total_energy(temperature, phonopy_dos[0], phonopy_dos[1])
+
+        harmonic_properties = [free_energy, entropy, c_v, total_energy, integration]
 
         # Renormalized force constants
         free_energy = thm.get_free_energy(temperature, phonopy_dos_r[0], phonopy_dos_r[1]) + \
-                      thm.get_free_energy_correction_2(temperature, phonopy_dos[0], phonopy_dos_r[1], phonopy_dos[1])
+                      thm.get_free_energy_correction_dos(temperature, phonopy_dos[0], phonopy_dos[1], phonopy_dos_r[1])
         entropy = thm.get_entropy(temperature, phonopy_dos_r[0], phonopy_dos_r[1])
         c_v = thm.get_cv(temperature, phonopy_dos_r[0], phonopy_dos_r[1])
+        total_energy = thm.get_total_energy(temperature, phonopy_dos_r[0], phonopy_dos_r[1]) + \
+                       thm.get_free_energy_correction_dos(temperature, phonopy_dos[0], phonopy_dos[1], phonopy_dos_r[1])
+
         integration = np.trapz(phonopy_dos_r[1], x=phonopy_dos_r[0])/(self.dynamic.structure.get_number_of_atoms()*
                                                        self.dynamic.structure.get_number_of_dimensions())
-        renormalized_properties = [free_energy, entropy, c_v, integration]
-        print('Free energy correction: {0:12.4f} KJ/mol'.format(thm.get_free_energy_correction_2(temperature, phonopy_dos[0], phonopy_dos_r[1], phonopy_dos[1])))
+        renormalized_properties = [free_energy, entropy, c_v, total_energy, integration]
+        print('Free energy correction: {0:12.4f} KJ/mol'.format(thm.get_free_energy_correction_dos(temperature, phonopy_dos[0], phonopy_dos[1], phonopy_dos_r[1])))
 
         if from_power_spectrum:
             normalization = np.prod(self.dynamic.get_super_cell_matrix())
 
-            power_spectrum = thm.get_dos(temperature, frequency_range, self.get_power_spectrum_direct(), normalization)
-            free_energy = thm.get_free_energy(temperature, frequency_range, power_spectrum)
-            entropy = thm.get_entropy(temperature, frequency_range, power_spectrum)
-            c_v = thm.get_cv(temperature, frequency_range, power_spectrum)
+            power_spectrum_dos = thm.get_dos(temperature, frequency_range, self.get_power_spectrum_direct(), normalization)
+            free_energy = thm.get_free_energy(temperature, frequency_range, power_spectrum_dos)
+            entropy = thm.get_entropy(temperature, frequency_range, power_spectrum_dos)
+            c_v = thm.get_cv(temperature, frequency_range, power_spectrum_dos)
+            total_energy = thm.get_total_energy(temperature, frequency_range, power_spectrum_dos)
 
-            integration = np.trapz(power_spectrum, x=frequency_range)/(self.dynamic.structure.get_number_of_atoms()*
+            integration = np.trapz(power_spectrum_dos, x=frequency_range)/(self.dynamic.structure.get_number_of_atoms()*
                                                            self.dynamic.structure.get_number_of_dimensions())
 
-            power_spectrum_properties = [free_energy, entropy, c_v, integration]
+            power_spectrum_properties = [free_energy, entropy, c_v, total_energy, integration]
             print('\nThermal properties per unit cell ({0:.2f} K) [From DoS]\n'
               '----------------------------------------------').format(temperature)
-            print('                         Harmonic    Renormalized   Power spectrum\n')
-            print('Free energy (KJ/mol): {0:12.4f}  {4:12.4f}  {8:12.4f}\n'
-                  'Entropy    (J/K/mol): {1:12.4f}  {5:12.4f}  {9:12.4f}\n'
-                  'Cv         (J/K/mol): {2:12.4f}  {6:12.4f}  {10:12.4f}\n'
-                  'Integration:          {3:12.4f}  {7:12.4f}  {11:12.4f}\n'.format(*(harmonic_properties +
+            print('                             Harmonic    Renormalized   Power spectrum\n')
+            print('Free energy   (KJ/mol): {0:12.4f}  {5:12.4f}  {10:12.4f}\n'
+                  'Entropy      (J/K/mol): {1:12.4f}  {6:12.4f}  {11:12.4f}\n'
+                  'Cv           (J/K/mol): {2:12.4f}  {7:12.4f}  {12:12.4f}\n'
+                  'Total energy  (KJ/mol): {3:12.4f}  {8:12.4f}  {13:12.4f}\n'
+                  'Integration:            {4:12.4f}  {9:12.4f}  {14:12.4f}\n'.format(*(harmonic_properties +
                                                                                       renormalized_properties +
                                                                                       power_spectrum_properties)))
             if not self.parameters.silent:
-                plt.plot(frequency_range, power_spectrum, 'r-', label='Molecular dynamics')
+                plt.plot(frequency_range, power_spectrum_dos, 'r-', label='Molecular dynamics')
 
         else:
             print('\nThermal properties per unit cell ({0:.2f} K) [From DoS]\n'
               '----------------------------------------------').format(temperature)
-            print('                        Harmonic    Renormalized\n')
-            print('Free energy  (KJ/mol): {0:12.4f}  {4:12.4f}\n'
-                  'Entropy     (J/K/mol): {1:12.4f}  {5:12.4f}\n'
-                  'Cv          (J/K/mol): {2:12.4f}  {6:12.4f}\n'
-                  'Integration:           {3:12.4f}  {7:12.4f}\n'.format(*(harmonic_properties + renormalized_properties)))
+            print('                            Harmonic    Renormalized\n')
+            print('Free energy   (KJ/mol): {0:12.4f}  {5:12.4f}\n'
+                  'Entropy      (J/K/mol): {1:12.4f}  {6:12.4f}\n'
+                  'Cv           (J/K/mol): {2:12.4f}  {7:12.4f}\n'
+                  'Total energy  (KJ/mol): {3:12.4f}  {8:12.4f}\n'
+                  'Integration:            {4:12.4f}  {9:12.4f}\n'.format(*(harmonic_properties + renormalized_properties)))
+
 
         if not self.parameters.silent:
             plt.plot(phonopy_dos[0], phonopy_dos[1], 'b-',label='Harmonic aprox.')

dynaphopy/analysis/thermal_properties.py

@@ -26,6 +26,18 @@ def energy2(freq, temp):
                          for i, freq in enumerate(frequency)])
     return dos
 
+def get_total_energy(temperature, frequency, dos):
+
+    def n(temp, freq):
+        return pow(np.exp(freq*h_bar/(k_b*temp))-1, -1)
+
+    total_energy = np.nan_to_num([dos[i] * h_bar*freq*(0.5 + n(temperature, freq))
+                                 for i, freq in enumerate(frequency)])
+
+    total_energy = np.trapz(total_energy, frequency) * N_a / 1000 # KJ/K/mol
+    return total_energy
+
+
 def get_free_energy(temperature, frequency, dos):
 
     free_energy = np.nan_to_num([dos[i] * k_b * temperature * np.log(2 * np.sinh(h_bar * freq / (2 * k_b * temperature)))
@@ -34,7 +46,7 @@ def get_free_energy(temperature, frequency, dos):
     free_energy = np.trapz(free_energy, frequency) * N_a / 1000 # KJ/K/mol
     return free_energy
 
-def get_free_energy_correction(temperature, frequency, dos, shift):
+def get_free_energy_correction_shift(temperature, frequency, dos, shift):
 
     def n(temp, freq):
         return pow(np.exp(freq*h_bar/(k_b*temp))-1, -1)
@@ -46,7 +58,7 @@ def n(temp, freq):
     return free_energy_c
 
 
-def get_free_energy_correction_2(temperature, frequency, dos, dos_r):
+def get_free_energy_correction_dos(temperature, frequency, dos, dos_r):
 
     def n(temp, freq):
         return pow(np.exp(freq*h_bar/(k_b*temp))-1, -1)
@@ -68,18 +80,20 @@ def get_entropy(temperature, frequency, dos):
     def coth(x):
         return  np.cosh(x)/np.sinh(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))))
+    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 = np.trapz(entropy, frequency) * N_a # J/K/mol
     return entropy
 
-#Alternative way to calculate entropy
+#Alternative way to calculate entropy (also works)
 def get_entropy2(temperature, frequency, dos):
 
     def n(temp, freq):
         return pow(np.exp(freq*h_bar/(k_b*temp))-1, -1)
 
-    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)))
+    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 = np.trapz(entropy, frequency) * N_a # J/K/mol
     return entropy
@@ -140,10 +154,10 @@ def z(temp, freq):
 
     #free_energy = get_free_energy(temp,frequency,dos) + get_free_energy_correction(temp, frequency, dos, shift)
 
-    print (get_free_energy_correction(temp, frequency, dos, shift),
-           get_free_energy_correction_2(temp, frequency, dos, dos_r))
+    print (get_free_energy_correction_shift(temp, frequency, dos, shift),
+           get_free_energy_correction_dos(temp, frequency, dos, dos_r))
 
-    free_energy = get_free_energy(temp, frequency_r, dos_r) + get_free_energy_correction_2(temp, frequency, dos_r, dos)
+    free_energy = get_free_energy(temp, frequency_r, dos_r) + get_free_energy_correction_dos(temp, frequency, dos_r, dos)
     entropy = get_entropy(temp, frequency_r, dos_r)
     c_v = get_cv(temp, frequency_r, dos_r)
     print ('Renormalized')

dynaphopy/classes/atoms.py

@@ -43,6 +43,7 @@ def __init__(self,
 
         self._super_cell_matrix = None
         self._super_cell_phonon = None
+        self._super_cell_phonon_renormalized = None
         self._number_of_cell_atoms = None
         self._number_of_atoms = None
         self._number_of_atom_types = None
@@ -101,9 +102,19 @@ def set_super_cell_phonon(self,super_cell_phonon):
 
     def get_super_cell_phonon(self):
         if self._super_cell_phonon is None:
-            self._super_cell_phonon = np.identity(self.get_number_of_dimensions(),dtype=int)
+            self._super_cell_phonon = np.identity(self.get_number_of_dimensions(), dtype=int)
         return self._super_cell_phonon
 
+    def set_super_cell_phonon_renormalized(self, super_cell_phonon):
+        self._super_cell_phonon_renormalized = super_cell_phonon
+
+
+    def get_super_cell_phonon_renormalized(self):
+        if self._super_cell_phonon_renormalized is None:
+            self._super_cell_phonon_renormalized = self.get_super_cell_phonon()
+        return self._super_cell_phonon_renormalized
+
+
 
     def set_super_cell_matrix(self,super_cell_matrix):
         self._super_cell_matrix = super_cell_matrix

dynaphopy/classes/parameters.py

@@ -17,7 +17,7 @@ def __init__(self,
 
                  # Fourier transform Method
                  correlation_function_step=10,
-                 integration_method = 1,  # 0: Trapezoid  1:Rectangles
+                 integration_method=1,  # 0: Trapezoid  1:Rectangles
 
                  # Fast Fourier tranform Method
                  zero_padding=0,
@@ -33,24 +33,27 @@ def __init__(self,
                  frequency_range=np.arange(0, 40.05, 0.05),
 
                  # Phonon dispersion diagram
-                 use_NAC = False,
+                 use_NAC=False,
                  band_ranges=([[[0.0, 0.0, 0.0], [0.5, 0.0, 0.5]]]),
-                 number_of_bins_histogram = 30,
+                 number_of_bins_histogram=30,
 
                  # Force constants
-                 symmetrize = False,
-                 use_symmetry = True,
+                 symmetrize=False,
+                 use_symmetry=True,
 
                  # Modes (eigenvectors) display
                  modes_vectors_scale=10,
 
                  #Density of states mesh (phonopy)
-                 mesh_phonopy=(40, 40, 40)
+                 mesh_phonopy=(40, 40, 40),
+
+                 #Use supercell
+                 use_MD_cell_commensurate=False,
                  ):
 
         self._silent = silent
-        self._number_of_coefficients_mem=number_of_coefficients_mem
-        self._mem_scan_range=mem_scan_range
+        self._number_of_coefficients_mem = number_of_coefficients_mem
+        self._mem_scan_range = mem_scan_range
         self._correlation_function_step = correlation_function_step
         self._integration_method = integration_method
         self._power_spectra_algorithm = power_spectra_algorithm
@@ -68,7 +71,7 @@ def __init__(self,
 
         self._modes_vectors_scale = modes_vectors_scale
         self._mesh_phonopy = mesh_phonopy
-
+        self._use_MD_cell_commensurate = use_MD_cell_commensurate
 
     def get_data_from_dict(self, data_dictionary):
         for data in self.__dict__:
@@ -222,3 +225,11 @@ def mesh_phonopy(self):
     @mesh_phonopy.setter
     def mesh_phonopy(self, mesh_phonopy):
         self._mesh_phonopy = mesh_phonopy
+
+    @property
+    def use_MD_cell_commensurate(self):
+        return self._use_MD_cell_commensurate
+
+    @use_MD_cell_commensurate.setter
+    def use_MD_cell_commensurate(self, use_MD_cell_commensurate):
+        self._use_MD_cell_commensurate = use_MD_cell_commensurate