Merge branch 'development'

Author: abelcarreras

README

@@ -1,5 +1,5 @@
 ---------------------------------------------------------
-DYNAPHOPY 1.11
+DYNAPHOPY 1.12
 ---------------------------------------------------------
 Software to calculate crystal microscopical anharmonic properties
 from molecular dynamics using the normal-mode-decomposition technique.

dynaphopy/__init__.py

@@ -1,4 +1,4 @@
-__version__='1.11'
+__version__='1.12'
 
 
 import numpy as np
@@ -15,6 +15,8 @@
 import dynaphopy.analysis.coordinates as trajdist
 import dynaphopy.analysis.thermal_properties as thm
 
+from scipy import integrate
+
 power_spectrum_functions = {
     0: [power_spectrum.get_fourier_spectra_par_openmp, 'Fourier transform'],
     1: [power_spectrum.get_mem_spectra_par_openmp,     'Maximum entropy method'],
@@ -282,7 +284,7 @@ def get_vc(self):
 
     def get_vq(self):
         if self._vq is None:
-            print("Projecting into phonon")
+            print("Projecting into phonon mode")
             self._vq =  projection.project_onto_phonon(self.get_vc(),self.get_eigenvectors())
         return self._vq
 
@@ -318,9 +320,9 @@ def plot_vc(self,atoms=None,coordinates=None):
         plt.legend()
         plt.show()
 
-    def save_vc(self,file_name):
+    def save_vc(self,file_name, atom=0):
         print("Saving wave vector projection to file")
-        np.savetxt(file_name, self.get_vc()[:, 0, :].real)
+        np.savetxt(file_name, self.get_vc()[:, atom, :].real)
 
     def save_vq(self,file_name):
         print("Saving phonon projection to file")
@@ -355,8 +357,8 @@ def get_power_spectrum_phonon(self):
                                                                                self.dynamic,
                                                                                self.parameters))
 
+                self.set_reduced_q_vector(initial_reduced_q_point)
                 self._power_spectrum_phonon = np.average(power_spectrum_phonon, axis=0)
-                self.parameters.reduced_q_vector = initial_reduced_q_point
             else:
                 self._power_spectrum_phonon = (
                     power_spectrum_functions[self.parameters.power_spectra_algorithm])[0](self.get_vq(),
@@ -383,8 +385,9 @@ def get_power_spectrum_wave_vector(self):
                                                                                            self.parameters))
 
                 power_spectrum_wave_vector = np.array(power_spectrum_wave_vector)
+                self.set_reduced_q_vector(initial_reduced_q_point)
                 self._power_spectrum_wave_vector = np.average(power_spectrum_wave_vector, axis=0)
-                self.parameters.reduced_q_vector = initial_reduced_q_point
+
             else:
                 self._power_spectrum_wave_vector = (
                         power_spectrum_functions[self.parameters.power_spectra_algorithm])[0](self.get_vc().swapaxes(1, 2).reshape(-1, size),
@@ -439,21 +442,26 @@ def phonon_individual_analysis(self):
 
     def plot_power_spectrum_full(self):
 
-        phonopy_dos = pho_interface.obtain_phonopy_dos(self.dynamic.structure,
-                                                       mesh=self.parameters.mesh_phonopy,
-                                                       freq_min=self.get_frequency_range()[0],
-                                                       freq_max=self.get_frequency_range()[-1],
-                                                       projected_on_atom=self.parameters.project_on_atom)
+        print(self.dynamic.structure.get_force_constants())
 
         fig, ax1 = plt.subplots()
 
         ax1.plot(self.get_frequency_range(), self.get_power_spectrum_direct(), 'r-', label='Power spectrum (MD)')
         ax1.set_xlabel('Frequency [THz]')
         ax1.set_ylabel('eV * ps')
-
         ax2 = ax1.twinx()
-        ax2.plot(phonopy_dos[0], phonopy_dos[1],'b-',label='DoS (Lattice dynamics)')
-        ax2.set_ylabel('Density of states')
+
+        if self.dynamic.structure.forces_available():
+            phonopy_dos = pho_interface.obtain_phonopy_dos(self.dynamic.structure,
+                                                           mesh=self.parameters.mesh_phonopy,
+                                                           freq_min=self.get_frequency_range()[0],
+                                                           freq_max=self.get_frequency_range()[-1],
+                                                           projected_on_atom=self.parameters.project_on_atom)
+
+
+
+            ax2.plot(phonopy_dos[0], phonopy_dos[1],'b-',label='DoS (Lattice dynamics)')
+            ax2.set_ylabel('Density of states')
 
         if self._renormalized_force_constants is not None:
             phonopy_dos_r = pho_interface.obtain_phonopy_dos(self.dynamic.structure,
@@ -475,7 +483,7 @@ def plot_power_spectrum_full(self):
 #        plt.legend()
         plt.show()
 
-        total_integral = np.trapz(self.get_power_spectrum_direct(), x=self.get_frequency_range())/(2 * np.pi)
+        total_integral = integrate.simps(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 plot_power_spectrum_wave_vector(self):
@@ -484,7 +492,7 @@ def plot_power_spectrum_wave_vector(self):
         plt.xlabel('Frequency [THz]')
         plt.ylabel('eV * ps')
         plt.show()
-        total_integral = np.trapz(self.get_power_spectrum_wave_vector(), x=self.get_frequency_range())/(2 * np.pi)
+        total_integral = integrate.simps(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))
 
 
@@ -571,14 +579,14 @@ 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)
+        total_integral = integrate.simps(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)
+        total_integral = integrate.simps(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):
@@ -755,7 +763,7 @@ def display_thermal_properties(self, from_power_spectrum=False, normalize_dos=Fa
         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 = np.trapz(phonopy_dos[1], x=phonopy_dos[0])/(self.dynamic.structure.get_number_of_atoms()*
+        integration = integrate.simps(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])
 
@@ -769,7 +777,7 @@ def display_thermal_properties(self, from_power_spectrum=False, normalize_dos=Fa
         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()*
+        integration = integrate.simps(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, total_energy, integration]
         print('Free energy/total 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])))
@@ -778,7 +786,7 @@ def display_thermal_properties(self, from_power_spectrum=False, normalize_dos=Fa
             normalization = np.prod(self.dynamic.get_super_cell_matrix())
 
             power_spectrum_dos = thm.get_dos(temperature, frequency_range, self.get_power_spectrum_direct(), normalization)
-            integration = np.trapz(power_spectrum_dos, x=frequency_range)/(self.dynamic.structure.get_number_of_atoms()*
+            integration = integrate.simps(power_spectrum_dos, x=frequency_range)/(self.dynamic.structure.get_number_of_atoms()*
                                                            self.dynamic.structure.get_number_of_dimensions())
 
             if normalize_dos:
@@ -823,7 +831,7 @@ def display_thermal_properties(self, from_power_spectrum=False, normalize_dos=Fa
             plt.plot(phonopy_dos[0], phonopy_dos[1], 'b-', label='Harmonic')
             plt.plot(phonopy_dos_r[0], phonopy_dos_r[1], 'g-', label='Renormalized')
 
-            plt.title('Density of states')
+            plt.title('Density of states (Normalized to unit cell)')
             plt.xlabel('Frequency [THz]')
             plt.ylabel('Density of states')
             plt.legend()

dynaphopy/analysis/coordinates.py

@@ -24,7 +24,8 @@ def progress_bar(progress):
 
 #print(disp.relative_trajectory(cell, traj ,pos))
 
-def relativize_trajectory(dynamic):
+#Not used (deprecated)
+def relativize_trajectory(dynamic, memmap=False):
 
     cell = dynamic.get_super_cell()
     number_of_atoms = dynamic.trajectory.shape[1]
@@ -35,13 +36,19 @@ def relativize_trajectory(dynamic):
 
     trajectory = dynamic.trajectory
 
+    if memmap:
+        normalized_trajectory = np.memmap('/home/abel/r_trajectory.map', dtype='complex', mode='w+', shape=trajectory.shape)
+    else:
+        normalized_trajectory = dynamic.trajectory.copy()
+
 #    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)
     return normalized_trajectory
+#Not used anymore (deprecated)
 
 def relativize_trajectory_py(dynamic):
 

dynaphopy/analysis/fitting.py

@@ -1,6 +1,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.optimize import curve_fit, minimize_scalar, root
+from scipy.integrate import simps
 
 h_planck = 4.135667662e-3  # eV/ps
 h_planck_bar = 6.58211951e-4  # eV/ps
@@ -129,7 +130,7 @@ def phonon_fitting_analysis(original, test_frequencies_range, harmonic_frequenci
         area = fit_params[2] / ( 2 * np.pi)
         base_line = fit_params[3]
 
-        total_integral = np.trapz(power_spectrum, x=test_frequencies_range)/ (2 * np.pi)
+        total_integral = simps(power_spectrum, x=test_frequencies_range)/ (2 * np.pi)
 
         #Calculated properties
         dt_Q2_lor = 2 * 2 * area

dynaphopy/analysis/thermal_properties.py

@@ -1,6 +1,7 @@
 import numpy as np
 import matplotlib.pylab as pl
 import warnings
+from scipy import integrate
 
 N_a = 6.022140857e23
 k_b = 1.38064852e-23  # J /K
@@ -36,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 = np.trapz(total_energy, frequency) * N_a / 1000  # KJ/K/mol
+    total_energy = integrate.simps(total_energy, frequency) * N_a / 1000  # KJ/K/mol
     return total_energy
 
 
@@ -46,7 +47,7 @@ def get_free_energy(temperature, frequency, dos):
                                  for i, freq in enumerate(frequency)])
 
     free_energy[0] = 0
-    free_energy = np.trapz(free_energy, frequency) * N_a / 1000  # KJ/K/mol
+    free_energy = integrate.simps(free_energy, frequency) * N_a / 1000  # KJ/K/mol
     return free_energy
 
 
@@ -58,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 = np.trapz(free_energy_c, frequency) * N_a / 1000 # KJ/K/mol
+    free_energy_c = integrate.simps(free_energy_c, frequency) * N_a / 1000 # KJ/K/mol
     return free_energy_c
 
 
@@ -75,7 +76,7 @@ def n(temp, freq):
 
     free_energy_c = free_energy_1 - free_energy_2
 
-    free_energy_c = np.trapz(free_energy_c, frequency) * N_a / 1000 # KJ/K/mol
+    free_energy_c = integrate.simps(free_energy_c, frequency) * N_a / 1000 # KJ/K/mol
     return free_energy_c
 
 
@@ -87,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 = np.trapz(entropy, frequency) * N_a # J/K/mol
+    entropy = integrate.simps(entropy, frequency) * N_a # J/K/mol
     return entropy
 
 #Alternative way to calculate entropy (also works)
@@ -99,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 = np.trapz(entropy, frequency) * N_a # J/K/mol
+    entropy = integrate.simps(entropy, frequency) * N_a # J/K/mol
     return entropy
 
 
@@ -110,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 = np.trapz(c_v, frequency) * N_a # J/K/mol
+    c_v = integrate.simps(c_v, frequency) * N_a # J/K/mol
 
     return c_v
 
@@ -170,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 ('\nFrom MD')
     print ('-------------------------')
@@ -182,6 +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 ('\nHARMONIC')
@@ -193,4 +195,5 @@ def z(temp, freq):
     print ('Free energy: {0} KJ/K/mol'.format(free_energy))
     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 (np.trapz(dos, x=frequency)/(8*3))
+    print (integrate.simps(dos, x=frequency)/(8*3))

dynaphopy/classes/atoms.py

@@ -189,9 +189,15 @@ def get_scaled_positions(self, supercell=None):
         return self._scaled_positions
 
 
+    #Force related methods
+
+    def forces_available(self):
+        if np.array(self.get_force_constants()).any() or np.array(self.get_force_sets()).any():
+            return True
+        else:
+            return False
 
 
-    #Force related methods
     def set_force_constants(self, force_constants):
         self._force_constants = np.array(force_constants,dtype='double')