Keep the consistency between the public ver. and the private ver.

anphon/anharmonic_core.cpp

@@ -358,8 +358,24 @@ void AnharmonicCore::prepare_group_of_force_constants(const std::vector<FcsArray
     }
 }
 
-
 std::complex<double> AnharmonicCore::V3(const unsigned int ks[3])
+{
+    return V3(ks,
+              dynamical->eval_phonon,
+              dynamical->evec_phonon);
+}
+
+std::complex<double> AnharmonicCore::V4(const unsigned int ks[4])
+{
+    return V4(ks,
+              dynamical->eval_phonon,
+              dynamical->evec_phonon);
+}
+
+
+std::complex<double> AnharmonicCore::V3(const unsigned int ks[3],
+                                        double **eval_phonon,
+                                        std::complex<double> ***evec_phonon)
 {
     int i;
     unsigned int kn[3], sn[3];
@@ -373,7 +389,7 @@ std::complex<double> AnharmonicCore::V3(const unsigned int ks[3])
     for (i = 0; i < 3; ++i) {
         kn[i] = ks[i] / ns;
         sn[i] = ks[i] % ns;
-        omega[i] = dynamical->eval_phonon[kn[i]][sn[i]];
+        omega[i] = eval_phonon[kn[i]][sn[i]];
     }
 
     // Return zero if any of the involving phonon has imaginary frequency
@@ -388,9 +404,9 @@ std::complex<double> AnharmonicCore::V3(const unsigned int ks[3])
 #pragma omp parallel for private(ret), reduction(+: ret_re, ret_im)
 #endif
     for (i = 0; i < ngroup_v3; ++i) {
-        ret = dynamical->evec_phonon[kn[0]][sn[0]][evec_index_v3[i][0]]
-            * dynamical->evec_phonon[kn[1]][sn[1]][evec_index_v3[i][1]]
-            * dynamical->evec_phonon[kn[2]][sn[2]][evec_index_v3[i][2]]
+        ret = evec_phonon[kn[0]][sn[0]][evec_index_v3[i][0]]
+            * evec_phonon[kn[1]][sn[1]][evec_index_v3[i][1]]
+            * evec_phonon[kn[2]][sn[2]][evec_index_v3[i][2]]
             * invmass_v3[i] * phi3_reciprocal[i];
         ret_re += ret.real();
         ret_im += ret.imag();
@@ -495,7 +511,9 @@ void AnharmonicCore::calc_phi3_reciprocal(const unsigned int ik1,
 }
 
 
-std::complex<double> AnharmonicCore::V4(const unsigned int ks[4])
+std::complex<double> AnharmonicCore::V4(const unsigned int ks[4],
+                                        double **eval_phonon,
+                                        std::complex<double> ***evec_phonon)
 {
     int i;
     int ns = dynamical->neval;
@@ -508,7 +526,7 @@ std::complex<double> AnharmonicCore::V4(const unsigned int ks[4])
     for (i = 0; i < 4; ++i) {
         kn[i] = ks[i] / ns;
         sn[i] = ks[i] % ns;
-        omega[i] = dynamical->eval_phonon[kn[i]][sn[i]];
+        omega[i] = eval_phonon[kn[i]][sn[i]];
     }
     // Return zero if any of the involving phonon has imaginary frequency
     if (omega[0] < eps8 || omega[1] < eps8 || omega[2] < eps8 || omega[3] < eps8) return 0.0;
@@ -531,10 +549,10 @@ std::complex<double> AnharmonicCore::V4(const unsigned int ks[4])
 #pragma omp parallel for private(ret), reduction(+: ret_re, ret_im)
 #endif
     for (i = 0; i < ngroup_v4; ++i) {
-        ret = dynamical->evec_phonon[kn[0]][sn[0]][evec_index_v4[i][0]]
-            * dynamical->evec_phonon[kn[1]][sn[1]][evec_index_v4[i][1]]
-            * dynamical->evec_phonon[kn[2]][sn[2]][evec_index_v4[i][2]]
-            * dynamical->evec_phonon[kn[3]][sn[3]][evec_index_v4[i][3]]
+        ret = evec_phonon[kn[0]][sn[0]][evec_index_v4[i][0]]
+            * evec_phonon[kn[1]][sn[1]][evec_index_v4[i][1]]
+            * evec_phonon[kn[2]][sn[2]][evec_index_v4[i][2]]
+            * evec_phonon[kn[3]][sn[3]][evec_index_v4[i][3]]
             * invmass_v4[i] * phi4_reciprocal[i];
         ret_re += ret.real();
         ret_im += ret.imag();
@@ -806,7 +824,9 @@ void AnharmonicCore::calc_damping_smearing(const unsigned int N,
             arr[1] = ns * k1 + is;
             arr[2] = ns * k2 + js;
 
-            v3_arr[ik][ib] = std::norm(V3(arr)) * multi;
+            v3_arr[ik][ib] = std::norm(V3(arr,
+                                          dynamical->eval_phonon,
+                                          dynamical->evec_phonon)) * multi;
         }
     }
 
@@ -1003,7 +1023,9 @@ void AnharmonicCore::calc_damping_tetrahedron(const unsigned int N,
                 arr[1] = ns * k1 + is;
                 arr[2] = ns * k2 + js;
 
-                v3_arr[ik][ib] = std::norm(V3(arr));
+                v3_arr[ik][ib] = std::norm(V3(arr,
+                                              dynamical->eval_phonon,
+                                              dynamical->evec_phonon));
 
             } else {
                 v3_arr[ik][ib] = 0.0;
@@ -1352,7 +1374,9 @@ void AnharmonicCore::calc_self3omega_tetrahedron(const double Temp,
                 arr[1] = ns * kpairs[ik_now][0] + is;
                 arr[2] = ns * kpairs[ik_now][1] + js;
 
-                v3_arr_loc[ib] = std::norm(V3(arr));
+                v3_arr_loc[ib] = std::norm(V3(arr,
+                                              dynamical->eval_phonon,
+                                              dynamical->evec_phonon));
             }
         }
         MPI_Gather(&v3_arr_loc[0], ns2, MPI_DOUBLE,

anphon/anharmonic_core.h

@@ -119,6 +119,14 @@ namespace PHON_NS
         std::complex<double> V3(const unsigned int [3]);
         std::complex<double> V4(const unsigned int [4]);
 
+        std::complex<double> V3(const unsigned int [3],
+                                double **,
+                                std::complex<double> ***);
+
+        std::complex<double> V4(const unsigned int [4],
+                                double **,
+                                std::complex<double> ***);
+
         std::complex<double> V3_mode(int,
                                      double *,
                                      double *,

anphon/conductivity.cpp

@@ -186,38 +186,40 @@ void Conductivity::prepare_restart()
         writes->fs_result.open(writes->file_result.c_str(), std::ios::app | std::ios::out);
     }
 
-
     // Add vks_done list here
 
     if (mympi->my_rank == 0) {
         nks_done = vks_done.size();
     }
     MPI_Bcast(&nks_done, 1, MPI_INT, 0, MPI_COMM_WORLD);
-    memory->allocate(arr_done, nks_done);
+    nshift_restart = nks_done;
 
-    if (mympi->my_rank == 0) {
-        for (i = 0; i < nks_done; ++i) {
-            arr_done[i] = vks_done[i];
-        }
-    }
-    MPI_Bcast(arr_done, nks_done, MPI_INT, 0, MPI_COMM_WORLD);
+    if (nks_done > 0) {
+        memory->allocate(arr_done, nks_done);
 
-    nshift_restart = nks_done;
+        if (mympi->my_rank == 0) {
+            for (i = 0; i < nks_done; ++i) {
+                arr_done[i] = vks_done[i];
+            }
+        }
+        MPI_Bcast(&arr_done[0], nks_done, MPI_INT, 0, MPI_COMM_WORLD);
 
-    // Remove vks_done elements from vks_job
+        // Remove vks_done elements from vks_job
 
-    for (i = 0; i < nks_done; ++i) {
+        for (i = 0; i < nks_done; ++i) {
 
-        it_set = vks_job.find(arr_done[i]);
+            it_set = vks_job.find(arr_done[i]);
 
-        if (it_set == vks_job.end()) {
-            error->exit("prepare_restart", "This cannot happen");
-        } else {
-            vks_job.erase(it_set);
+            if (it_set == vks_job.end()) {
+                std::cout << " rank = " << mympi->my_rank
+                    << " arr_done = " << arr_done[i] << std::endl;
+                error->exit("prepare_restart", "This cannot happen");
+            } else {
+                vks_job.erase(it_set);
+            }
         }
+        memory->deallocate(arr_done);
     }
-
-    memory->deallocate(arr_done);
     vks_done.clear();
 }
 

anphon/dynamical.cpp

@@ -779,7 +779,7 @@ void Dynamical::diagonalize_dynamical_all()
         std::cout << "done!" << std::endl;
     }
 
-    if (kpoint->kpoint_mode == 2) {
+    if (kpoint->kpoint_mode == 2 && phon->mode == "RTA") {
         detect_imaginary_branches(dynamical->eval_phonon);
     }
 }

anphon/fcs_phonon.cpp

@@ -72,6 +72,7 @@ void Fcs_phonon::setup(std::string mode)
 
     MPI_Bcast(&anharmonic_core->quartic_mode, 1, MPI_INT, 0, MPI_COMM_WORLD);
     MPI_Bcast(&gruneisen->print_gruneisen, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD);
+    MPI_Bcast(&thermodynamics->calc_FE_bubble, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD);
 
     if (mode == "PHONONS") {
         require_cubic = false;

anphon/kpoint.h

@@ -206,7 +206,7 @@ namespace PHON_NS
         void get_unique_triplet_k(const int,
                                   const bool,
                                   const bool,
-                                  std::vector<KsListGroup> &, 
+                                  std::vector<KsListGroup> &,
                                   const int sign = -1);
 
     private:

anphon/mode_analysis.cpp

@@ -268,7 +268,7 @@ void ModeAnalysis::run_mode_analysis()
         print_selfenergy(NT, T_arr);
 
         if (print_V3) print_V3_elements();
-//        if (print_V3) print_Phi3_elements();
+        //        if (print_V3) print_Phi3_elements();
 
         if (calc_fstate_omega) print_frequency_resolved_final_state(NT, T_arr);
 
@@ -1850,13 +1850,12 @@ void ModeAnalysis::calc_Phi3(const unsigned int knum,
             omega[2] = dynamical->eval_phonon[k2][js];
 
             ret[ik][ib] = anharmonic_core->V3(arr) * factor
-            * std::sqrt(omega[0] * omega[1] * omega[2]);
+                * std::sqrt(omega[0] * omega[1] * omega[2]);
         }
     }
 }
 
 
-
 void ModeAnalysis::print_spectral_function(const int NT,
                                            double *T_arr)
 {

anphon/phonons.cpp

@@ -208,6 +208,10 @@ void PHON::execute_phonons()
         gruneisen->calc_gruneisen();
     }
 
+    if (thermodynamics->calc_FE_bubble) {
+        thermodynamics->compute_free_energy_bubble();
+    }
+
     if (mympi->my_rank == 0) {
         writes->print_phonon_energy();
         writes->write_phonon_info();