Diegom spectral function multiprocessing

cellconstructor/Phonons.py

@@ -1254,15 +1254,15 @@ def save_phononpy(self, *args, **kwargs):
         warnings.warn("[DEPRECATION WARNING] save_phononpy is deprecated: use save_phonopy instead.")
         self.save_phonopy(*args, **kwargs)
 
-    def save_phonopy(self, path = ".", supercell_size = None, units_ev_ang2 = True, 
+    def save_phonopy(self, path = ".", supercell_size = None, units_ev_ang2 = True,
         write_poscar = True, write_unitcell = False):
         """
         EXPORT THE DYN IN THE PHONONPY FORMAT
         =====================================
 
         This tool export the dynamical matrix into the PHONONPY plain text format.
         If units_ev_ang2 is True (default) the dynamical matrix is written in eV/A^2
-        Otherwise we use the Ry/bohr^2, as the quantum espresso format. 
+        Otherwise we use the Ry/bohr^2, as the quantum espresso format.
         Please, remember this when using Phononpy for the conversion factors.
 
         It will create a file called FORCE_CONSTANTS, one called unitcell.in
@@ -1276,13 +1276,13 @@ def save_phonopy(self, path = ".", supercell_size = None, units_ev_ang2 = True,
                 The supercell that defines the dynamical matrix, note phononpy
                 works in the supercell. If none, it is inferred from the q points
             units_ev_ang2 : bool
-                If True (default) convert the units in eV / A^2 
+                If True (default) convert the units in eV / A^2
             write_poscar : bool
                 If True produce also the POSCAR file with the structure
                 for phonopy.
                 It requires ASE to be installed.
             write_unitcell : bool
-                If true, produce a unitcell.in for phonopy. 
+                If true, produce a unitcell.in for phonopy.
 
 
         """
@@ -3494,7 +3494,7 @@ def DiagonalizeSupercell(self, verbose = False, lo_to_split = None, return_qmode
         DYAGONALIZE THE DYNAMICAL MATRIX IN THE SUPERCELL
         =================================================
 
-        This method dyagonalizes the dynamical matrix in q space 
+        This method dyagonalizes the dynamical matrix in q space
         and returns the frequencies and the polarization vectors in the supercell,
         without having to generate the force constant in real space.
 
@@ -4450,7 +4450,7 @@ def get_dyn_from_ase_phonons(ase_ph, adjust_qstar = True):
 
 
 
-def compute_phonons_finite_displacements(structure, ase_calculator, epsilon = 0.05, 
+def compute_phonons_finite_displacements(structure, ase_calculator, epsilon = 0.05,
     supercell = (1,1,1), progress = -1, progress_bar = False,
     use_symmetries = True, debug=False):
     """
@@ -4588,11 +4588,11 @@ def compute_force(indices):
 
 
 
-def compute_phonons_finite_displacements_sym(structure, ase_calculator, epsilon=0.05, 
+def compute_phonons_finite_displacements_sym(structure, ase_calculator, epsilon=0.05,
                                              supercell=(1,1,1),
                                              progress=-1,
                                              progress_bar=False,
-                                             debug=False, 
+                                             debug=False,
                                              timer=None):
     """
     COMPUTE THE FORCE CONSTANT MATRIX EXPLOITING SYMMETRIES
@@ -4668,7 +4668,7 @@ def compute_phonons_finite_displacements_sym(structure, ase_calculator, epsilon=
     # Get irt from the symmetries (atomic corrispondance after the application of symmetry)
     if debug:
         print("Getting symmetry equivalent atoms")
-    
+
     if timer is not None:
         irts = timer.execute_timed_function(symmetries.get_symmetry_equivalent_atoms, symm, super_structure)
     else:
@@ -4749,7 +4749,7 @@ def compute_force(indices):
 
     # Define the force constant matrix in the basis of the auxiliary vectors
     fc_aux_basis = np.zeros((nat3, nat3), dtype = np.double)
-    
+
     # Compute the auxiliary force basis
     # This could exploit parallelization to speedup the calculation
     for index, gen in enumerate(generators):
@@ -4766,7 +4766,7 @@ def compute_force(indices):
             force_sym = timer.execute_timed_function(symmetries.ApplySymmetryToVector,current_sym, force, super_structure.unit_cell, current_irt)
         else:
             force_sym = symmetries.ApplySymmetryToVector(current_sym, force, super_structure.unit_cell, current_irt)
-        
+
         fc_aux_basis[index, :] = force_sym.ravel()
 
 
@@ -4802,7 +4802,7 @@ def compute_force(indices):
 
     #                 if coeffs is None:
     #                     disp_basis.append(v.ravel())
-    #                     counter_index += 1                
+    #                     counter_index += 1
     #                     fc_aux_basis[counter_index, :] = force_sym[i_sym, :, :].ravel()
 
     #np.savetxt("OriginalFC.dat", fc, fmt="%10.6f")
@@ -4835,7 +4835,7 @@ def compute_force(indices):
         correct_dyn = Phonons(structure, nqirr = np.prod(supercell))
         q_tot = symmetries.GetQGrid(structure.unit_cell, supercell)
         if timer is not None:
-            dynq = timer.execute_timed_function(GetDynQFromFCSupercell, 
+            dynq = timer.execute_timed_function(GetDynQFromFCSupercell,
                 final_dyn.dynmats[0], np.array(q_tot), structure, super_structure)
         else:
             dynq = GetDynQFromFCSupercell(final_dyn.dynmats[0], np.array(q_tot), structure, super_structure)

cellconstructor/Spectral.py

@@ -241,11 +241,10 @@ def get_static_correction_interpolated(dyn, tensor3, T, new_supercell, k_grid):
     dynq = get_static_correction(dyn, tensor3, k_grid, q_tot, T)
 
     # Add all the new computed dynamical matrix
-    all_dynmats = []
-    for iq in range(len(q_tot)):
-        all_dynmats.append(dynq[iq, :, :])
-
-    new_dyn.dynmats = all_dynmats
+    new_dyn.dynmats[0] = dynq[0, :, :]
+    for iq in range(1,len(q_tot)):
+        new_dyn.dynmats.append(dynq[iq, :, :])
+
 
     # Adjust the dynamical matrix q points and the stars
     new_dyn.AdjustQStar()
@@ -1178,20 +1177,22 @@ def output_file_sort_function(filename_sp, smear_id_cm, smear_cm, nsm):
             data = np.loadtxt(f)
             f.close()
             head = head1[1:]+head2[1:]+head3[1:]
-            # X = data[:,0]
-            # Y = data[:,1]
-            # Z = data[:,2]
-            # x = [X[i] for i in np.lexsort((Y,X))]
-            # y = [Y[i] for i in np.lexsort((Y,X))]
-            # z = [Z[i] for i in np.lexsort((Y,X))]
+            X = data[:,0]
+            Y = data[:,1]
+            Z = data[:,2]
+            x = [X[i] for i in np.lexsort((Y,X))]
+            y = [Y[i] for i in np.lexsort((Y,X))]
+            z = [Z[i] for i in np.lexsort((Y,X))]
             # f = open('plot_energies_'+filename_data, 'w')
             # np.savetxt(f,np.c_[x,y,z])
             # f.close()
-            data = data[data[:,2].argsort()] # First sort doesn't need to be stable.
-            data = data[data[:,1].argsort(kind='mergesort')]
-            data = data[data[:,0].argsort(kind='mergesort')]
+            # data = data[data[:,2].argsort()] # First sort doesn't need to be stable.
+            # data = data[data[:,1].argsort(kind='mergesort')]
+            # data = data[data[:,0].argsort(kind='mergesort')]
             f = open("Sorted_"+filename_data, 'w')
-            np.savetxt(f,data, header=head[:-1])
+            #np.savetxt(f,data, header=head[:-1])
+            np.savetxt(f,np.c_[x,y,z], header=head[:-1])
+            f.close()
         pass
 
     print(" ")
@@ -1331,16 +1332,20 @@ def output_file_sort_function(filename_sp, smear_id_cm, smear_cm, nsm):
     plwork.join()
 
 
-
     if static_limit :
         print(" ")
         print(" Results printed in "+filename_sp+'_static.dat')
         print(" ")
+        filename_new=filename_sp+'_static.dat'
+        output_file_sort_function(filename_new, smear_id_cm, smear_cm, nsm)
     else:
         print(" ")
         print(" Results printed in "+filename_sp+'_[smear].dat')
         print(" ")
-    output_file_sort_function(filename_sp, smear_id_cm, smear_cm, nsm)
+        for  ism in range(nsm):
+            name="{:6.1f}".format(smear_cm[ism]).strip()
+            filename_new=filename_sp+'_'+name+'.dat'
+            output_file_sort_function(filename_sp, smear_id_cm, smear_cm, nsm)
 
 
 
@@ -3472,7 +3477,7 @@ def get_os_perturb_dynamic_correction_along_path(dyn, tensor3,
         q_path : list of triplets
                  Path of the q-points of the Brillouin Zone, in 2pi/Anstrom units,
                  where the caculation is performed
-                 (defualt: [0.0,0.0,0.0])
+                 (default: [0.0,0.0,0.0])
         q_path_file : string
                       Name of the file where the q_path can be read.
                       Format: column of triples, q points in 2pi/Angstrom