1919#include " mpi_common.h"
2020#include " phonon_dos.h"
2121#include " pointers.h"
22+ #include " progress_bar.h"
2223#include " relaxation.h"
2324#include " system.h"
24- #include " progress_bar.h"
2525
2626using namespace PHON_NS ;
2727
28- Thermodynamics::Thermodynamics (PHON *phon) : Pointers(phon)
29- {
30- T_to_Ryd = k_Boltzmann / Ryd;
31- calc_FE_bubble = false ;
32- FE_bubble = nullptr ;
33- }
28+ Thermodynamics::Thermodynamics (PHON *phon) :
29+ Pointers(phon), T_to_Ryd(k_Boltzmann / Ryd), classical(false ), calc_FE_bubble(false ), FE_bubble(nullptr )
30+ {}
3431
3532Thermodynamics::~Thermodynamics ()
3633{
@@ -44,40 +41,40 @@ void Thermodynamics::setup()
4441 MPI_Bcast (&classical, 1 , MPI_CXX_BOOL, 0 , MPI_COMM_WORLD);
4542}
4643
47- double Thermodynamics::Cv (const double omega, const double temp_in) const
44+ auto Thermodynamics::Cv (const double omega, const double temp_in) const -> double
4845{
4946 if (std::abs (temp_in) < eps) return 0.0 ;
5047
5148 const auto x = omega / (T_to_Ryd * temp_in);
5249 return k_Boltzmann * std::pow (x / (2.0 * sinh (0.5 * x)), 2 );
5350}
5451
55- double Thermodynamics::Cv_classical (const double omega, const double temp_in) const
52+ auto Thermodynamics::Cv_classical (const double omega, const double temp_in) -> double
5653{
5754 if (std::abs (temp_in) < eps) return 0.0 ;
5855
5956 return k_Boltzmann;
6057}
6158
62- double Thermodynamics::fB (const double omega, const double temp_in) const
59+ auto Thermodynamics::fB (const double omega, const double temp_in) const -> double
6360{
6461 if (std::abs (temp_in) < eps || omega < eps8) return 0.0 ;
6562
6663 const auto x = omega / (T_to_Ryd * temp_in);
6764 return 1.0 / (std::exp (x) - 1.0 );
6865}
6966
70- double Thermodynamics::fC (const double omega, const double temp_in) const
67+ auto Thermodynamics::fC (const double omega, const double temp_in) const -> double
7168{
7269 if (std::abs (temp_in) < eps || omega < eps8) return 0.0 ;
7370
7471 const auto x = omega / (T_to_Ryd * temp_in);
7572 return 1.0 / x;
7673}
7774
78- double Thermodynamics::Cv_tot (const double temp_in, const unsigned int nk_irred, const unsigned int ns,
79- const std::vector<std::vector<KpointList>> &kp_irred, const double *weight_k_irred,
80- const double *const *eval_in) const
75+ auto Thermodynamics::Cv_tot (const double temp_in, const unsigned int nk_irred, const unsigned int ns,
76+ const std::vector<std::vector<KpointList>> &kp_irred, const double *weight_k_irred,
77+ const double *const *eval_in) const -> double
8178{
8279 int i;
8380 unsigned int ik, is;
@@ -118,10 +115,10 @@ double Thermodynamics::Cv_tot(const double temp_in, const unsigned int nk_irred,
118115 return ret;
119116}
120117
121- double Thermodynamics::Cv_anharm_correction (const double temp_in, const unsigned int nk_irred, const unsigned int ns,
122- const std::vector<std::vector<KpointList>> &kp_irred,
123- const double *weight_k_irred, const double *const *eval_in,
124- const double *const *del_eval_in) const
118+ auto Thermodynamics::Cv_anharm_correction (const double temp_in, const unsigned int nk_irred, const unsigned int ns,
119+ const std::vector<std::vector<KpointList>> &kp_irred,
120+ const double *weight_k_irred, const double *const *eval_in,
121+ const double *const *del_eval_in) const -> double
125122{
126123 int i;
127124 unsigned int ik, is;
@@ -164,9 +161,9 @@ double Thermodynamics::Cv_anharm_correction(const double temp_in, const unsigned
164161 return ret;
165162}
166163
167- double Thermodynamics::internal_energy (const double temp_in, const unsigned int nk_irred, const unsigned int ns,
168- const std::vector<std::vector<KpointList>> &kp_irred,
169- const double *weight_k_irred, const double * const * eval_in) const
164+ auto Thermodynamics::internal_energy (const double temp_in, const unsigned int nk_irred, const unsigned int ns,
165+ const std::vector<std::vector<KpointList>> &kp_irred, const double *weight_k_irred ,
166+ const double *const * eval_in) const -> double
170167{
171168 int i;
172169 unsigned int ik, is;
@@ -206,9 +203,9 @@ double Thermodynamics::internal_energy(const double temp_in, const unsigned int
206203 return ret * 0.5 ;
207204}
208205
209- double Thermodynamics::vibrational_entropy (const double temp_in, const unsigned int nk_irred, const unsigned int ns,
210- const std::vector<std::vector<KpointList>> &kp_irred,
211- const double *weight_k_irred, const double *const *eval_in) const
206+ auto Thermodynamics::vibrational_entropy (const double temp_in, const unsigned int nk_irred, const unsigned int ns,
207+ const std::vector<std::vector<KpointList>> &kp_irred,
208+ const double *weight_k_irred, const double *const *eval_in) const -> double
212209{
213210 int i;
214211 unsigned int ik, is;
@@ -249,9 +246,9 @@ double Thermodynamics::vibrational_entropy(const double temp_in, const unsigned
249246 return -k_Boltzmann * ret;
250247}
251248
252- double Thermodynamics::free_energy_QHA (const double temp_in, const unsigned int nk_irred, const unsigned int ns,
253- const std::vector<std::vector<KpointList>> &kp_irred,
254- const double *weight_k_irred, const double * const * eval_in) const
249+ auto Thermodynamics::free_energy_QHA (const double temp_in, const unsigned int nk_irred, const unsigned int ns,
250+ const std::vector<std::vector<KpointList>> &kp_irred, const double *weight_k_irred ,
251+ const double *const * eval_in) const -> double
255252{
256253 int i;
257254 unsigned int ik, is;
@@ -302,25 +299,25 @@ double Thermodynamics::free_energy_QHA(const double temp_in, const unsigned int
302299 return temp_in * T_to_Ryd * ret;
303300}
304301
305- double Thermodynamics::disp2_avg (const double T_in, const unsigned int ncrd1, const unsigned int ncrd2,
306- const unsigned int nk, const unsigned int ns, const double *const *xk_in,
307- const double *const *eval_in, std::complex <double > ***evec_in) const
302+ auto Thermodynamics::disp2_avg (const double T_in, const unsigned int ncrd1, const unsigned int ncrd2,
303+ const unsigned int nk, const unsigned int ns, const double *const *xk_in,
304+ const double *const *eval_in, std::complex <double > ***evec_in) const -> double
308305{
309- const double cell_shift[3 ]{0 , 0 , 0 };
306+ constexpr double cell_shift[3 ]{0 , 0 , 0 };
310307 return disp_corrfunc (T_in, ncrd1, ncrd2, cell_shift, nk, ns, xk_in, eval_in, evec_in);
311308}
312309
313- double Thermodynamics::disp_corrfunc (const double T_in, const unsigned int ncrd1, const unsigned int ncrd2,
314- const double cell_shift[3 ], const unsigned int nk, const unsigned int ns,
315- const double *const *xk_in, const double *const *eval_in,
316- std::complex <double > ***evec_in) const
310+ auto Thermodynamics::disp_corrfunc (const double T_in, const unsigned int ncrd1, const unsigned int ncrd2,
311+ const double cell_shift[3 ], const unsigned int nk, const unsigned int ns,
312+ const double *const *xk_in, const double *const *eval_in,
313+ std::complex <double > ***evec_in) const -> double
317314{
318315 int i;
319- int N = nk * ns;
316+ int const N = nk * ns;
320317 unsigned int ik, is;
321318 double omega;
322319 double phase;
323- const std::complex <double > im (0.0 , 1.0 );
320+ constexpr std::complex <double > im (0.0 , 1.0 );
324321 double ret = 0.0 ;
325322
326323 if (classical) {
@@ -362,7 +359,7 @@ double Thermodynamics::disp_corrfunc(const double T_in, const unsigned int ncrd1
362359 return ret;
363360}
364361
365- double Thermodynamics::coth_T (const double omega, const double T) const
362+ auto Thermodynamics::coth_T (const double omega, const double T) const -> double
366363{
367364 // This function returns coth(hbar*omega/2*kB*T)
368365
@@ -373,7 +370,7 @@ double Thermodynamics::coth_T(const double omega, const double T) const
373370 return 1.0 + 2.0 / (std::exp (2.0 * x) - 1.0 );
374371}
375372
376- void Thermodynamics::compute_free_energy_bubble ()
373+ auto Thermodynamics::compute_free_energy_bubble () -> void
377374{
378375 const auto NT = static_cast <unsigned int >((system->Tmax - system->Tmin ) / system->dT ) + 1 ;
379376
@@ -392,7 +389,7 @@ void Thermodynamics::compute_free_energy_bubble()
392389 }
393390}
394391
395- void Thermodynamics::compute_FE_bubble (double **eval, std::complex <double > ***evec, double *FE_bubble_out) const
392+ auto Thermodynamics::compute_FE_bubble (double **eval, std::complex <double > ***evec, double *FE_bubble_out) const -> void
396393{
397394 // This function calculates the free energy of the bubble diagram
398395 double omega_sum[2 ];
@@ -451,7 +448,7 @@ void Thermodynamics::compute_FE_bubble(double **eval, std::complex<double> ***ev
451448 triplet,
452449 1 );
453450
454- const int npair_uniq = triplet.size ();
451+ const size_t npair_uniq = triplet.size ();
455452
456453 arr_cubic[0 ] = ns * ik0 + is0;
457454
@@ -523,24 +520,21 @@ void Thermodynamics::compute_FE_bubble(double **eval, std::complex<double> ***ev
523520 deallocate (FE_tmp);
524521}
525522
526- void Thermodynamics::compute_FE_bubble_SCPH (double ***eval_in, std::complex <double > ****evec_in, double *FE_bubble)
523+ auto Thermodynamics::compute_FE_bubble_SCPH (double ***eval_in, std::complex <double > ****evec_in,
524+ double *FE_bubble) const -> void
527525{
528526 // This function calculates the free energy from the bubble diagram
529527 // at the given temperature and lattice dynamics wavefunction
530- int ik;
531- int multi;
532- double omega0, omega1, omega2, omega_sum[2 ];
528+ double omega_sum[2 ];
533529 double n0, n1, n2, nsum[2 ];
534530 const auto nk = dos->kmesh_dos ->nk ;
535531 const auto nk_reduced = dos->kmesh_dos ->nk_irred ;
536532 const auto ns = dynamical->neval ;
537- double v3_tmp;
538- unsigned int ik0, ik1, ik2, is0, is1, is2, i0, iT;
533+ unsigned int i0, iT;
539534 unsigned int arr_cubic[3 ];
540535 const auto nks0 = nk_reduced * ns;
541- unsigned int NT = static_cast <unsigned int >((system->Tmax - system->Tmin ) / system->dT ) + 1 ;
542- double temp;
543- double factor = -1.0 / (static_cast <double >(nk * nk) * 48.0 );
536+ const unsigned int NT = static_cast <unsigned int >((system->Tmax - system->Tmin ) / system->dT ) + 1 ;
537+ const double factor = -1.0 / (static_cast <double >(nk * nk) * 48.0 );
544538
545539 double *FE_local;
546540 double *FE_tmp;
@@ -558,9 +552,9 @@ void Thermodynamics::compute_FE_bubble_SCPH(double ***eval_in, std::complex<doub
558552 }
559553 }
560554
561- auto startTime = std::chrono::system_clock::now ();
555+ const auto startTime = std::chrono::system_clock::now ();
562556 auto lastUpdate = startTime;
563- bool isConsole = isOutputToConsole ();
557+ const bool isConsole = isOutputToConsole ();
564558
565559 unsigned int nk_tmp;
566560
@@ -573,11 +567,11 @@ void Thermodynamics::compute_FE_bubble_SCPH(double ***eval_in, std::complex<doub
573567 vks_l.push_back (-1 );
574568 }
575569
576- auto nks_tmp = vks_l.size ();
570+ const auto nks_tmp = vks_l.size ();
577571
578-
579- for (iT = 0 ; iT < NT; ++iT)
572+ for (iT = 0 ; iT < NT; ++iT) {
580573 FE_local[iT] = 0.0 ;
574+ }
581575
582576 if (mympi->my_rank == 0 ) {
583577 std::cout << " Total number of modes per MPI process: " << nk_tmp << ' \n ' ;
@@ -587,8 +581,8 @@ void Thermodynamics::compute_FE_bubble_SCPH(double ***eval_in, std::complex<doub
587581
588582 if (vks_l[i0] != -1 ) {
589583
590- ik0 = dos->kmesh_dos ->kpoint_irred_all [vks_l[i0] / ns][0 ].knum ;
591- is0 = vks_l[i0] % ns;
584+ unsigned int const ik0 = dos->kmesh_dos ->kpoint_irred_all [vks_l[i0] / ns][0 ].knum ;
585+ unsigned int const is0 = vks_l[i0] % ns;
592586
593587 dos->kmesh_dos ->get_unique_triplet_k (vks_l[i0] / ns,
594588 symmetry->SymmList ,
@@ -597,43 +591,44 @@ void Thermodynamics::compute_FE_bubble_SCPH(double ***eval_in, std::complex<doub
597591 triplet,
598592 1 );
599593
600- int npair_uniq = triplet.size ();
594+ const size_t npair_uniq = triplet.size ();
601595
602596 arr_cubic[0 ] = ns * ik0 + is0;
603597
604598 for (iT = 0 ; iT < NT; ++iT)
605599 FE_tmp[iT] = 0.0 ;
606600
607- for (ik = 0 ; ik < npair_uniq; ++ik) {
608- multi = static_cast <double >(triplet[ik].group .size ());
601+ for (size_t ik = 0 ; ik < npair_uniq; ++ik) {
602+ const auto multi = static_cast <double >(triplet[ik].group .size ());
609603
610604 arr_cubic[0 ] = ns * ik0 + is0;
611605
612- ik1 = triplet[ik].group [0 ].ks [0 ];
613- ik2 = triplet[ik].group [0 ].ks [1 ];
606+ const unsigned int ik1 = triplet[ik].group [0 ].ks [0 ];
607+ const unsigned int ik2 = triplet[ik].group [0 ].ks [1 ];
614608
615- for (is1 = 0 ; is1 < ns; ++is1) {
609+ for (unsigned int is1 = 0 ; is1 < ns; ++is1) {
616610 arr_cubic[1 ] = ns * ik1 + is1;
617611
618- for (is2 = 0 ; is2 < ns; ++is2) {
612+ for (unsigned int is2 = 0 ; is2 < ns; ++is2) {
619613 arr_cubic[2 ] = ns * ik2 + is2;
620614
621615 for (iT = 0 ; iT < NT; ++iT) {
622616
623- temp = system->Tmin + static_cast <double >(iT) * system->dT ;
617+ const double temp = system->Tmin + static_cast <double >(iT) * system->dT ;
624618
625- omega0 = eval_in[iT][ik0][is0];
626- omega1 = eval_in[iT][ik1][is1];
627- omega2 = eval_in[iT][ik2][is2];
619+ const double omega0 = eval_in[iT][ik0][is0];
620+ const double omega1 = eval_in[iT][ik1][is1];
621+ const double omega2 = eval_in[iT][ik2][is2];
628622
629623 if (omega0 < eps8 || omega1 < eps8 || omega2 < eps8) continue ;
630624
631625 omega_sum[0 ] = 1.0 / (omega0 + omega1 + omega2);
632626 omega_sum[1 ] = 1.0 / (-omega0 + omega1 + omega2);
633627
634- v3_tmp = std::norm (
635- anharmonic_core->V3 (arr_cubic, dos->kmesh_dos ->xk , eval_in[iT], evec_in[iT])) *
636- static_cast <double >(multi);
628+ const double v3_tmp =
629+ std::norm (
630+ anharmonic_core->V3 (arr_cubic, dos->kmesh_dos ->xk , eval_in[iT], evec_in[iT])) *
631+ multi;
637632
638633 if (classical) {
639634 n0 = fC (omega0, temp);
@@ -656,16 +651,17 @@ void Thermodynamics::compute_FE_bubble_SCPH(double ***eval_in, std::complex<doub
656651 }
657652 }
658653 }
659- double weight = static_cast <double >(dos->kmesh_dos ->kpoint_irred_all [vks_l[i0] / ns].size ());
660- for (iT = 0 ; iT < NT; ++iT)
654+ const double weight = static_cast <double >(dos->kmesh_dos ->kpoint_irred_all [vks_l[i0] / ns].size ());
655+ for (iT = 0 ; iT < NT; ++iT) {
661656 FE_local[iT] += FE_tmp[iT] * weight;
657+ }
662658 }
663659 if (mympi->my_rank == 0 ) {
664660 auto currentTime = std::chrono::system_clock::now ();
665- long long totalElapsedTime =
661+ const long long totalElapsedTime =
666662 std::chrono::duration_cast<std::chrono::milliseconds>(currentTime - startTime).count ();
667- long long avgTimePerStep = (i0 == 0 ) ? 0 : totalElapsedTime / i0;
668- long long timeRemaining = (i0 == 0 ) ? 0 : avgTimePerStep * (nks_tmp - i0 - 1 );
663+ const long long avgTimePerStep = (i0 == 0 ) ? 0 : totalElapsedTime / i0;
664+ const long long timeRemaining = (i0 == 0 ) ? 0 : avgTimePerStep * (nks_tmp - i0 - 1 );
669665 displayProgressBar (i0, nks_tmp - 1 , std::cout, timeRemaining, isConsole, " Fe-bubble" );
670666 lastUpdate = currentTime;
671667 if (i0 == nk_tmp - 1 ) std::cout << " \n done. \n\n " << std::flush;
@@ -682,9 +678,9 @@ void Thermodynamics::compute_FE_bubble_SCPH(double ***eval_in, std::complex<doub
682678 deallocate (FE_tmp);
683679}
684680
685- double Thermodynamics::FE_scph_correction (unsigned int iT, double **eval, std::complex <double > ***evec,
686- double **eval_harm_renormalized,
687- std::complex <double > ***evec_harm_renormalized) const
681+ auto Thermodynamics::FE_scph_correction (unsigned int iT, double **eval, std::complex <double > ***evec,
682+ double **eval_harm_renormalized,
683+ std::complex <double > ***evec_harm_renormalized) const -> double
688684{
689685 using namespace Eigen ;
690686 const auto nk = dos->kmesh_dos ->nk ;
@@ -735,7 +731,8 @@ double Thermodynamics::FE_scph_correction(unsigned int iT, double **eval, std::c
735731 return ret / static_cast <double >(nk);
736732}
737733
738- double Thermodynamics::compute_FE_total (unsigned int iT, double fe_qha, double dfe_scph = 0.0 )
734+ auto Thermodynamics::compute_FE_total (const unsigned int iT, const double fe_qha, const double dfe_scph = 0.0 ) const
735+ -> double
739736{
740737 double fe_total = fe_qha;
741738 // skip scph correction for QHA + structural optimization
0 commit comments