print out the quasiparticle entropy for numerical evaluation of Cv

ttadano · ttadano · commit 80b09073cf93 · 2026-01-16T17:20:31.000+09:00
diff --git a/anphon/scph.cpp b/anphon/scph.cpp
@@ -295,6 +295,7 @@ void Scph::postprocess(std::complex<double> ****delta_dymat, std::complex<double
         double *FE_QHA = nullptr;
         double *dFE_scph = nullptr;
         double *FE_total = nullptr;
+        double *entropy = nullptr;
         double **msd_update = nullptr;
         double ***ucorr_update = nullptr;
         double ****dielec_update = nullptr;
@@ -318,6 +319,7 @@ void Scph::postprocess(std::complex<double> ****delta_dymat, std::complex<double
             allocate(FE_QHA, NT);
             allocate(dFE_scph, NT);
             allocate(FE_total, NT);
+            allocate(entropy, NT);
 
             if (writes->getPrintMSD()) {
                 allocate(msd_update, NT, ns);
@@ -434,6 +436,14 @@ void Scph::postprocess(std::complex<double> ****delta_dymat, std::complex<double
 
                 FE_total[iT] = thermodynamics->compute_FE_total(iT, FE_QHA[iT], dFE_scph[iT]);
 
+                entropy[iT] = thermodynamics->vibrational_entropy(temperature,
+                                                                  dos->kmesh_dos->nk_irred,
+                                                                  ns,
+                                                                  dos->kmesh_dos->kpoint_irred_all,
+                                                                  dos->kmesh_dos->weight_k.data(),
+                                                                  eval_update[iT]) /
+                              k_Boltzmann;
+
                 if (writes->getPrintMSD()) {
                     double shift[3]{0.0, 0.0, 0.0};
 
@@ -501,7 +511,13 @@ void Scph::postprocess(std::complex<double> ****delta_dymat, std::complex<double
             if (dos->compute_dos) {
                 writes->writePhononDos(dos_update, is_qha, 0);
             }
-            writes->writeThermodynamicFunc(heat_capacity, heat_capacity_correction, FE_QHA, dFE_scph, FE_total, is_qha);
+            writes->writeThermodynamicFunc(heat_capacity,
+                                           heat_capacity_correction,
+                                           FE_QHA,
+                                           dFE_scph,
+                                           FE_total,
+                                           entropy,
+                                           is_qha);
             if (writes->getPrintMSD()) {
                 writes->writeMSD(msd_update, is_qha, 0);
             }
@@ -777,6 +793,7 @@ void Scph::postprocess(std::complex<double> ****delta_dymat, std::complex<double
         if (FE_QHA) deallocate(FE_QHA);
         if (dFE_scph) deallocate(dFE_scph);
         if (FE_total) deallocate(FE_total);
+        if (entropy) deallocate(entropy);
         if (dielec_update) deallocate(dielec_update);
 
         if (eval_gam) deallocate(eval_gam);
@@ -2131,8 +2148,6 @@ void Scph::interpolate_to_dense_mesh(std::complex<double> ***dymat_q,
 }
 
 
-
-
 bool Scph::check_convergence(const Eigen::MatrixXd &omega_now, const Eigen::MatrixXd &omega_old, const double conv_tol,
                              const unsigned int verbosity, const int iloop, double &diff) const
 {
@@ -2386,10 +2401,10 @@ void Scph::compute_anharmonic_frequency(std::complex<double> ***v4_array_all, do
 
 
 void Scph::compute_anharmonic_frequency_diis_perkpoint(std::complex<double> ***v4_array_all, double **omega2_out,
-                                                        std::complex<double> ***evec_anharm_scph, const double temp,
-                                                        bool &flag_converged, std::complex<double> ***cmat_convert,
-                                                        const bool offdiag, std::complex<double> **delta_v2_renorm,
-                                                        const unsigned int verbosity)
+                                                       std::complex<double> ***evec_anharm_scph, const double temp,
+                                                       bool &flag_converged, std::complex<double> ***cmat_convert,
+                                                       const bool offdiag, std::complex<double> **delta_v2_renorm,
+                                                       const unsigned int verbosity)
 {
     // SCPH with per-k-point DIIS using GDIIS_PerKpoint class
     // Each k-point has independent DIIS history and coefficients
@@ -2448,13 +2463,23 @@ void Scph::compute_anharmonic_frequency_diis_perkpoint(std::complex<double> ***v
 
     const auto T_in = temp;
 
-    initialize_scph_iteration(T_in, flag_converged, omega2_out, verbosity,
-                              omega_now, omega2_HA, evec_initial, evec_initial_adjoint, cmat_convert);
+    initialize_scph_iteration(T_in,
+                              flag_converged,
+                              omega2_out,
+                              verbosity,
+                              omega_now,
+                              omega2_HA,
+                              evec_initial,
+                              evec_initial_adjoint,
+                              cmat_convert);
 
     setup_harmonic_dynamical_matrices(omega2_HA, evec_initial, delta_v2_renorm, Fmat0, dymat_q_HA);
 
-    dynamical->precompute_dymat_harm(kmesh_dense->nk, kmesh_dense->xk, kmesh_dense->kvec_na,
-                                     dymat_harm_short, dymat_harm_long);
+    dynamical->precompute_dymat_harm(kmesh_dense->nk,
+                                     kmesh_dense->xk,
+                                     kmesh_dense->kvec_na,
+                                     dymat_harm_short,
+                                     dymat_harm_long);
 
     // Initialize per-k-point DIIS mixer
     const int diis_history = std::min(3, static_cast<int>(maxiter / 3));
@@ -2473,13 +2498,28 @@ void Scph::compute_anharmonic_frequency_diis_perkpoint(std::complex<double> ***v
             knum = kmap_interpolate_to_scph[knum_interpolate];
 
             update_fmat_with_v4(Fmat0, v4_array_all, dmat_convert, offdiag, ik, Fmat);
-            diagonalize_and_symmetrize(Fmat, evec_initial, v4_array_all, ik, knum, knum_interpolate,
-                                       flag_converged, omega2_out, verbosity, icount, eval_tmp, dymat_q);
+            diagonalize_and_symmetrize(Fmat,
+                                       evec_initial,
+                                       v4_array_all,
+                                       ik,
+                                       knum,
+                                       knum_interpolate,
+                                       flag_converged,
+                                       omega2_out,
+                                       verbosity,
+                                       icount,
+                                       eval_tmp,
+                                       dymat_q);
         }
 
         // Step 2: Interpolate to get new eigenvalues (output from using dmat_convert)
-        interpolate_to_dense_mesh(dymat_q, dymat_q_HA, evec_initial, eval_interpolate,
-                                  evec_new, cmat_convert, omega_now);
+        interpolate_to_dense_mesh(dymat_q,
+                                  dymat_q_HA,
+                                  evec_initial,
+                                  eval_interpolate,
+                                  evec_new,
+                                  cmat_convert,
+                                  omega_now);
 
         // Step 3: Check convergence before mixing
         if (check_convergence(omega_now, omega_old, conv_tol, verbosity, iloop, diff)) {
@@ -2492,7 +2532,7 @@ void Scph::compute_anharmonic_frequency_diis_perkpoint(std::complex<double> ***v
         }
 
         if (use_diis) {
-            gdiis_mixer_perkpoint.push(dmat_convert, omega_now);  // Correct: D_input → ω_output
+            gdiis_mixer_perkpoint.push(dmat_convert, omega_now); // Correct: D_input → ω_output
         }
 
         // Step 5: Compute what D should be based on new eigenvalues
@@ -2549,7 +2589,8 @@ void Scph::compute_anharmonic_frequency_diis_perkpoint(std::complex<double> ***v
                 std::cout << " (with per-k-point DIIS)";
                 auto stats = gdiis_mixer_perkpoint.get_success_stats();
                 int total_success = 0;
-                for (int s : stats) total_success += s;
+                for (int s: stats)
+                    total_success += s;
                 if (verbosity > 1) {
                     std::cout << "\n  Total DIIS successes: " << total_success;
                 }
diff --git a/anphon/thermodynamics.cpp b/anphon/thermodynamics.cpp
@@ -43,6 +43,9 @@ void Thermodynamics::setup()
 
 auto Thermodynamics::Cv(const double omega, const double temp_in) const -> double
 {
+    // Mode specific heat at constant volume.
+    // Exactly heat capacity only within the QHA.
+    // In other cases, Σ_q C_v is not equal to the total heat capacity.
     if (std::abs(temp_in) < eps) return 0.0;
 
     const auto x = omega / (T_to_Ryd * temp_in);
@@ -51,13 +54,15 @@ auto Thermodynamics::Cv(const double omega, const double temp_in) const -> doubl
 
 auto Thermodynamics::Cv_classical(const double omega, const double temp_in) -> double
 {
+    // Just return k_B
     if (std::abs(temp_in) < eps) return 0.0;
 
     return k_Boltzmann;
 }
 
 auto Thermodynamics::fB(const double omega, const double temp_in) const -> double
 {
+    // Bose-Einstein distribution function
     if (std::abs(temp_in) < eps || omega < eps8) return 0.0;
 
     const auto x = omega / (T_to_Ryd * temp_in);
@@ -66,6 +71,7 @@ auto Thermodynamics::fB(const double omega, const double temp_in) const -> doubl
 
 auto Thermodynamics::fC(const double omega, const double temp_in) const -> double
 {
+    // Classical limit of Bose-Einstein distribution function
     if (std::abs(temp_in) < eps || omega < eps8) return 0.0;
 
     const auto x = omega / (T_to_Ryd * temp_in);
@@ -76,6 +82,7 @@ auto Thermodynamics::Cv_tot(const double temp_in, const unsigned int nk_irred, c
                             const std::vector<std::vector<KpointList>> &kp_irred, const double *weight_k_irred,
                             const double *const *eval_in) const -> double
 {
+    // Total constant-volume heat capacity. Only the quasiharmonic term is included here.
     int i;
     unsigned int ik, is;
     double omega;
@@ -120,6 +127,10 @@ auto Thermodynamics::Cv_anharm_correction(const double temp_in, const unsigned i
                                           const double *weight_k_irred, const double *const *eval_in,
                                           const double *const *del_eval_in) const -> double
 {
+    // Anharmonic correction to constant-volume heat capacity
+    // We only use the adjacent temperature point for the numerical derivative, so the numerical accuracy
+    // may not be very high. For more reliable estimate, it is recommended to fit the entropy curve with
+    // respect to temperature by polynomial function and take the derivative of the fitted curve.
     int i;
     unsigned int ik, is;
     double omega, domega_dt;
@@ -165,6 +176,8 @@ auto Thermodynamics::internal_energy(const double temp_in, const unsigned int nk
                                      const std::vector<std::vector<KpointList>> &kp_irred, const double *weight_k_irred,
                                      const double *const *eval_in) const -> double
 {
+    // Vibrational internal energy within QHA.
+    // U = F + TS = Σ_q  0.5 * hw_q * coth(hw_q/2k_BT)] = Σ_q [hw_q (n_q + 0.5)]
     int i;
     unsigned int ik, is;
     double omega;
@@ -207,6 +220,8 @@ auto Thermodynamics::vibrational_entropy(const double temp_in, const unsigned in
                                          const std::vector<std::vector<KpointList>> &kp_irred,
                                          const double *weight_k_irred, const double *const *eval_in) const -> double
 {
+    // Vibrational entropy correct for QHA/SCP and other quasiparticle approaches
+    // S = -(∂F/∂T) = k_B * Σ_q [(n_q + 1)ln(n_q + 1) - n_q ln n_q]
     int i;
     unsigned int ik, is;
     double omega, x;
@@ -250,6 +265,8 @@ auto Thermodynamics::free_energy_QHA(const double temp_in, const unsigned int nk
                                      const std::vector<std::vector<KpointList>> &kp_irred, const double *weight_k_irred,
                                      const double *const *eval_in) const -> double
 {
+    // Vibrational free energy within QHA and QHA-like term within SCP.
+    // F = Σ_q [0.5 hw_q + k_B T ln(1 - exp(-hw_q/k_BT))] (quantum)
     int i;
     unsigned int ik, is;
     double omega, x;
@@ -682,6 +699,8 @@ auto Thermodynamics::FE_scph_correction(unsigned int iT, double **eval, std::com
                                         double **eval_harm_renormalized,
                                         std::complex<double> ***evec_harm_renormalized) const -> double
 {
+    // The correction term to the free energy within SCPH theory.
+    // This term is necessary to result in S =
     using namespace Eigen;
     const auto nk = dos->kmesh_dos->nk;
     const auto ns = dynamical->neval;
diff --git a/anphon/write_phonons.cpp b/anphon/write_phonons.cpp
@@ -3131,7 +3131,7 @@ void Writes::writePhononDos(double **dos_in, const bool is_qha, const int bubble
 }
 
 void Writes::writeThermodynamicFunc(double *heat_capacity, double *heat_capacity_correction, double *FE_QHA,
-                                    double *dFE_scph, double *FE_total, const bool is_qha) const
+                                    double *dFE_scph, double *FE_total, double *entropy, const bool is_qha) const
 {
     const auto Tmin = system->Tmin;
     const auto Tmax = system->Tmax;
@@ -3161,6 +3161,12 @@ void Writes::writeThermodynamicFunc(double *heat_capacity, double *heat_capacity
     }
     if (thermodynamics->calc_FE_bubble) {
         ofs_thermo << "# The bubble free-energy calculated on top of the SCPH wavefunction is also shown.\n";
+        ofs_thermo << "# However, the bubble contributions to the heat capacity and entropy are not included.\n";
+        ofs_thermo << "# If these are needed, please fit the free energy data including the bubble term \n "
+                      "# by polynomial function and then estimate S and Cv by numerical derivatives.\n";
+    }
+    if (!is_qha) {
+        ofs_thermo << "# The Cv data accounts for the QHA-like term only.\n";
     }
 
     ofs_thermo << "# Temperature [K], Cv [in kB unit]";
@@ -3175,21 +3181,11 @@ void Writes::writeThermodynamicFunc(double *heat_capacity, double *heat_capacity
     if (thermodynamics->calc_FE_bubble) {
         ofs_thermo << ", F_{vib} (Bubble correction) [Ry]";
     }
-    // write renormalized zero-th order IFC
+        // write renormalized zero-th order IFC
     if (relaxation->relax_str != 0) {
         ofs_thermo << ", Phi0 [Ry]";
     }
-    ofs_thermo << ", F_{total} [Ry]\n";
-
-    // if (thermodynamics->calc_FE_bubble) {
-    //     ofs_thermo << "# The bubble free-energy calculated on top of the SCPH wavefunction is also shown." << '\n';
-    //     ofs_thermo <<
-    //                "# Temperature [K], Cv [in kB unit], F_{vib} (QHA term) [Ry], F_{vib} (SCPH correction) [Ry], F_{vib} (Bubble correction) [Ry]"
-    //                << '\n';
-    // } else {
-    //     ofs_thermo << "# Temperature [K], Cv [in kB unit], F_{vib} (QHA term) [Ry], F_{vib} (SCPH correction) [Ry]"
-    //                << '\n';
-    // }
+    ofs_thermo << ", F_{total} [Ry], S_{vib} [in kB unit]\n";
 
     if (thermodynamics->classical) {
         ofs_thermo << "# CLASSICAL = 1: Use classical limit.\n";
@@ -3216,15 +3212,16 @@ void Writes::writeThermodynamicFunc(double *heat_capacity, double *heat_capacity
         if (relaxation->relax_str != 0) {
             ofs_thermo << std::setw(18) << relaxation->V0[iT];
         }
-        ofs_thermo << std::setw(18) << FE_total[iT] << '\n';
+        ofs_thermo << std::setw(18) << FE_total[iT];
+        ofs_thermo << std::setw(18) << entropy[iT] << '\n';
     }
 
     ofs_thermo.close();
     std::cout << "  " << std::setw(input->job_title.length() + 12) << std::left << file_thermo;
     if (is_qha) {
-        std::cout << " : QHA heat capcaity and free energy\n";
+        std::cout << " : QHA heat capcaity, free energy, entropy\n";
     } else {
-        std::cout << " : SCPH heat capcaity and free energy\n";
+        std::cout << " : SCPH heat capcaity, free energy, entropy\n";
     }
 }
 
diff --git a/anphon/write_phonons.h b/anphon/write_phonons.h
@@ -58,7 +58,7 @@ class Writes: protected Pointers
     void writePhononDos(double **dos_in, const bool is_qha = false, const int bubble = 0) const;
 
     void writeThermodynamicFunc(double *heat_capacity, double *heat_capacity_correction, double *FE_QHA,
-                                double *dFE_scph, double *FE_total, const bool is_qha = false) const;
+                                double *dFE_scph, double *FE_total, double *entropy, const bool is_qha = false) const;
 
     void writeMSD(double **msd_in, const bool is_qha = false, const int bubble = 0) const;
 
diff --git a/tools/plotband.py b/tools/plotband.py
@@ -43,10 +43,10 @@ def get_kpath_and_kval(file_in):
         kval_float = [float(val) for val in kval[1:]]
         kpath_list = []
         for i in range(len(kpath[1:])):
-            if kpath[i + 1] == "G":
-                kpath_list.append("$\Gamma$")
+            if kpath[i + 1] == "G" or kpath[i + 1] == "GAMMA":
+                kpath_list.append(r"$\Gamma$")
             else:
-                kpath_list.append("$\mathrm{%s}$" % kpath[i + 1])
+                kpath_list.append(r"$\mathrm{%s}$" % kpath[i + 1])
 
         return kpath_list, kval_float
     else:
