Implement ISOTOPE = 2 to save phonon selfenergy due to phonon-isotope…

… scatterings
ttadano · Jan 15, 2015 · 7817f16 · 7817f16
1 parent 611ecb5
commit 7817f16
Show file tree

Hide file tree

Showing 7 changed files with 62 additions and 5 deletions.
diff --git a/anphon/conductivity.cpp b/anphon/conductivity.cpp
@@ -409,6 +409,8 @@ void Conductivity::compute_kappa()
             }
         }
         memory->deallocate(lifetime);
+
+        std::cout << "" << std::endl;
     }
 }
 

diff --git a/anphon/isotope.cpp b/anphon/isotope.cpp
@@ -39,7 +39,7 @@ void Isotope::setup_isotope_scattering()
     int i;
     int nkd = system->nkd;
 
-    MPI_Bcast(&include_isotope, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD);
+    MPI_Bcast(&include_isotope, 1, MPI_INT, 0, MPI_COMM_WORLD);
 
     if (include_isotope) {
 
@@ -51,7 +51,7 @@ void Isotope::setup_isotope_scattering()
 
         if (mympi->my_rank == 0) {
             std::cout << std::endl;
-            std::cout << " ISOTOPE = 1: Isotope scattering effects will be considered" << std::endl;
+            std::cout << " ISOTOPE >= 1: Isotope scattering effects will be considered" << std::endl;
             std::cout << "              with the following scattering factors." << std::endl;
 
             for (i = 0; i < nkd; ++i) {
@@ -247,5 +247,7 @@ void Isotope::calc_isotope_selfenergy_all()
 
 */
     }
+
+
 }
 
diff --git a/anphon/isotope.h b/anphon/isotope.h
@@ -20,7 +20,7 @@ namespace PHON_NS
         Isotope(class PHON *);
         ~Isotope();
 
-        bool include_isotope;
+        int include_isotope;
         double *isotope_factor;
 
         double **gamma_isotope;

diff --git a/anphon/parsephon.cpp b/anphon/parsephon.cpp
@@ -286,7 +286,6 @@ void Input::parse_analysis_vars(const bool use_default_values)
                                    ISOTOPE ISOFACT FSTATE_W FSTATE_K PRINTMSD PDOS TDOS GRUNEISEN NEWFCS DELTA_A \
                                    ANIME ANIME_CELLSIZE ANIME_FORMAT SPS";
 
-    bool include_isotope;
     bool fstate_omega, fstate_k;
     bool lclassical;
     bool ks_analyze_mode, atom_project_mode, calc_realpart;
@@ -296,6 +295,8 @@ void Input::parse_analysis_vars(const bool use_default_values)
     bool print_xsf, print_anime;
 
     int quartic_mode;
+    int include_isotope;
+
 
     double delta_a;
     double *isotope_factor;
@@ -325,7 +326,7 @@ void Input::parse_analysis_vars(const bool use_default_values)
     ks_analyze_mode = false;
     atom_project_mode = false;
     calc_realpart = false;
-    include_isotope = false;
+    include_isotope = 0;
     fstate_omega = false;
     fstate_k = false;
 

diff --git a/anphon/phonons.cpp b/anphon/phonons.cpp
@@ -247,6 +247,7 @@ void PHON::execute_RTA()
         conductivity->calc_anharmonic_imagself();
         conductivity->compute_kappa();
         writes->write_kappa();
+        writes->write_selfenergy_isotope();
     }
 
     if (kpoint->kpoint_mode == 2) {

diff --git a/anphon/write_phonons.cpp b/anphon/write_phonons.cpp
@@ -1120,10 +1120,60 @@ void Writes::write_kappa()
         ofs_kl.close();
 
         std::cout << std::endl;
+        std::cout << " ------------------------------------------------------------" << std::endl << std::endl;
         std::cout << " Lattice thermal conductivity is store in the file " << file_kappa << std::endl;
     }
 }
 
+void Writes::write_selfenergy_isotope()
+{
+    unsigned int i, j, k;
+    unsigned int ns = dynamical->neval;
+    unsigned int knum;
+    double **eval = dynamical->eval_phonon;
+    double **gamma_iso = isotope->gamma_isotope;
+
+    if (mympi->my_rank == 0) {
+        if (isotope->include_isotope == 2) {
+
+            std::string file_iso = input->job_title + ".self_isotope";
+            std::ofstream ofs_iso;
+
+            ofs_iso.open(file_iso.c_str(), std::ios::out);
+            if (!ofs_iso) error->exit("write_selfenergy_isotope", "Could not open file_iso");
+
+            ofs_iso << "# Phonon selfenergy due to phonon-isotope scatterings for the irreducible k points." << std::endl;
+            ofs_iso << "# Irred. knum, mode num, frequency [cm^-1], Gamma_iso [cm^-1]" << std::endl << std::endl;
+
+            for (i = 0; i < kpoint->nk_reduced; ++i) {
+                ofs_iso << "# Irreducible k point  : " << std::setw(8) << i + 1;
+                ofs_iso << " (" << std::setw(4) << kpoint->kpoint_irred_all[i].size() << ")" << std::endl;
+
+                knum = kpoint->kpoint_irred_all[i][0].knum;
+
+                ofs_iso << "## xk = " << std::setw(3);
+                for (k = 0; k < 3; ++k) ofs_iso << std::setw(15) << kpoint->xk[knum][k];
+                ofs_iso << std::endl;
+
+                for (k = 0; k < ns; ++k){
+                    ofs_iso << std::setw(7) << i + 1;
+                    ofs_iso << std::setw(5) << k + 1;
+                    ofs_iso << std::setw(15) << in_kayser(eval[knum][k]);
+                    ofs_iso << std::setw(15) << in_kayser(gamma_iso[i][k]);
+                    ofs_iso << std::endl;
+                }
+                ofs_iso << std::endl;
+            }
+
+
+            std::cout << std::endl;
+            std::cout << " ISOTOPE = 2: Phonon selfenergy due to phonon-isotope " << std::endl;
+            std::cout << "              scatterings is stored in the file " << file_iso << std::endl;
+
+            ofs_iso.close();
+        }
+    }
+}
 
 void Writes::write_normal_mode_animation(const double xk_in[3], const unsigned int ncell[3])
 {

diff --git a/anphon/write_phonons.h b/anphon/write_phonons.h
@@ -29,6 +29,7 @@ namespace PHON_NS
         void setup_result_io();
         void write_input_vars();
         void write_kappa();
+        void write_selfenergy_isotope();
 
         bool print_xsf;
         bool print_msd;