Skip to content

Commit

Permalink
Implement classical version of FE_bubble and FE_SCPH
Browse files Browse the repository at this point in the history
  • Loading branch information
@ttadano
ttadano committed Nov 28, 2018
1 parent aefc874 commit cfadd78
Show file tree
Hide file tree
Showing 2 changed files with 76 additions and 24 deletions.
59 changes: 46 additions & 13 deletions anphon/scph.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -311,7 +311,7 @@ void Scph::load_scph_dymat_from_file(std::complex<double> ****dymat_out)
"The number of KMESH_SCPH is not consistent");
}
if (nonanalytic_tmp != dynamical->nonanalytic) {
error->exit("load_scph_dymat_from_file",
error->warn("load_scph_dymat_from_file",
"The NONANALYTIC tag is not consistent");
}
if (consider_offdiag_tmp != consider_offdiagonal) {
Expand Down Expand Up @@ -2954,42 +2954,70 @@ void Scph::write_scph_thermodynamics(double ***eval,
<< std::endl;
}

if (thermodynamics->classical) {
ofs_thermo << "# CLASSICAL = 1: Use classical limit." << std::endl;
}

for (unsigned int iT = 0; iT < NT; ++iT) {

double temp = Tmin + static_cast<double>(iT) * dT;

double tmp1 = 0.0;
double tmp2 = 0.0;
double heat_capacity = 0.0;
double free_energy = 0.0;

#pragma omp parallel for reduction(+:tmp1,tmp2)
if (thermodynamics->classical) {
#pragma omp parallel for reduction(+:heat_capacity,free_energy)
for (int i = 0; i < N; ++i) {
unsigned int ik = i / ns;
unsigned int is = i % ns;
double omega = eval[iT][ik][is];

if (omega <= eps8) continue;

tmp1 += thermodynamics->Cv(omega, temp);
heat_capacity += thermodynamics->Cv_classical(omega, temp);

if (std::abs(temp) > eps) {
double x = omega / (temp * T_to_Ryd);
free_energy += std::log(x);
}
}

heat_capacity /= static_cast<double>(nk);
free_energy *= temp * T_to_Ryd / static_cast<double>(nk);

} else {

#pragma omp parallel for reduction(+:heat_capacity,free_energy)
for (int i = 0; i < N; ++i) {
unsigned int ik = i / ns;
unsigned int is = i % ns;
double omega = eval[iT][ik][is];

if (omega <= eps8) continue;

heat_capacity += thermodynamics->Cv(omega, temp);

if (std::abs(temp) < eps) {
tmp2 += 0.5 * omega;
free_energy += 0.5 * omega;
} else {
double x = omega / (temp * T_to_Ryd);
tmp2 += 0.5 * x + std::log(1.0 - std::exp(-x));
free_energy += 0.5 * x + std::log(1.0 - std::exp(-x));
}
}

tmp1 /= static_cast<double>(nk);
heat_capacity /= static_cast<double>(nk);
if (std::abs(temp) < eps) {
tmp2 /= static_cast<double>(nk);
free_energy /= static_cast<double>(nk);
} else {
tmp2 *= temp * T_to_Ryd / static_cast<double>(nk);
free_energy *= temp * T_to_Ryd / static_cast<double>(nk);
}
double tmp3 = tmp2 + FE_scph(iT, eval[iT], evec[iT]);
}

double tmp3 = free_energy + FE_scph(iT, eval[iT], evec[iT]);

ofs_thermo << std::setw(16) << std::fixed << temp;
ofs_thermo << std::setw(18) << std::scientific << tmp1 / k_Boltzmann;
ofs_thermo << std::setw(18) << tmp2;
ofs_thermo << std::setw(18) << std::scientific << heat_capacity / k_Boltzmann;
ofs_thermo << std::setw(18) << free_energy;
ofs_thermo << std::setw(18) << tmp3;

if (thermodynamics->calc_FE_bubble) {
Expand Down Expand Up @@ -3038,8 +3066,13 @@ double Scph::FE_scph(unsigned int iT,
}
}

if (thermodynamics->classical) {
ret += (tmp_c.real() - omega * omega) * thermodynamics->fC(omega, temp) / (4.0 * omega);

} else {
ret += (tmp_c.real() - omega * omega) * (1.0 + 2.0 * thermodynamics->fB(omega, temp)) / (8.0 * omega);
}
}

return ret / static_cast<double>(nk);
}
Expand Down
41 changes: 30 additions & 11 deletions anphon/thermodynamics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -435,12 +435,22 @@ void Thermodynamics::compute_FE_bubble(double **eval,

for (iT = 0; iT < NT; ++iT) {
double temp = system->Tmin + static_cast<double>(iT) * system->dT;
double n0 = thermodynamics->fB(omega0, temp);
double n1 = thermodynamics->fB(omega1, temp);
double n2 = thermodynamics->fB(omega2, temp);

nsum[0] = (1.0 + n0) * (1.0 + n1 + n2) + n1 * n2;
nsum[1] = n0 * n1 - n1 * n2 + n2 * n0 + n0;
if (classical) {
double n0 = fC(omega0, temp);
double n1 = fC(omega1, temp);
double n2 = fC(omega2, temp);

nsum[0] = n0 * (n1 + n2) + n1 * n2;
nsum[1] = n0 * (n1 + n2) - n1 * n2;
} else {
double n0 = fB(omega0, temp);
double n1 = fB(omega1, temp);
double n2 = fB(omega2, temp);

nsum[0] = (1.0 + n0) * (1.0 + n1 + n2) + n1 * n2;
nsum[1] = n0 * n1 - n1 * n2 + n2 * n0 + n0;
}

FE_tmp[iT] += v3_tmp * (nsum[0] * omega_sum[0] + 3.0 * nsum[1] * omega_sum[1]);
}
Expand Down Expand Up @@ -556,12 +566,21 @@ double Thermodynamics::compute_FE_bubble_SCPH(const double temp,
double v3_tmp = std::norm(anharmonic_core->V3(arr_cubic, eval, evec))
* static_cast<double>(multi);

double n0 = thermodynamics->fB(omega0, temp);
double n1 = thermodynamics->fB(omega1, temp);
double n2 = thermodynamics->fB(omega2, temp);

nsum[0] = (1.0 + n0) * (1.0 + n1 + n2) + n1 * n2;
nsum[1] = n0 * n1 - n1 * n2 + n2 * n0 + n0;
if (classical) {
double n0 = fC(omega0, temp);
double n1 = fC(omega1, temp);
double n2 = fC(omega2, temp);

nsum[0] = n0 * (n1 + n2) + n1 * n2;
nsum[1] = n0 * (n1 + n2) - n1 * n2;
} else {
double n0 = fB(omega0, temp);
double n1 = fB(omega1, temp);
double n2 = fB(omega2, temp);

nsum[0] = (1.0 + n0) * (1.0 + n1 + n2) + n1 * n2;
nsum[1] = n0 * n1 - n1 * n2 + n2 * n0 + n0;
}

FE_tmp += v3_tmp * (nsum[0] * omega_sum[0] + 3.0 * nsum[1] * omega_sum[1]);
}
Expand Down

0 comments on commit cfadd78

Please sign in to comment.