Skip to content

Commit

Permalink
Performance improvement of the thermodynamic function calculation
Browse files Browse the repository at this point in the history
  • Loading branch information
@ttadano
ttadano committed Mar 2, 2018
1 parent 62ee00b commit 45a6c63
Show file tree
Hide file tree
Showing 2 changed files with 223 additions and 57 deletions.
279 changes: 222 additions & 57 deletions anphon/thermodynamics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,11 +22,10 @@ using namespace PHON_NS;
Thermodynamics::Thermodynamics(PHON *phon): Pointers(phon)
{
T_to_Ryd = k_Boltzmann / Ryd;
sum_over_irreducible_points = true;
}

Thermodynamics::~Thermodynamics()
{
};
Thermodynamics::~Thermodynamics() {};

void Thermodynamics::setup()
{
Expand Down Expand Up @@ -78,29 +77,67 @@ double Thermodynamics::Cv_tot(const double T)
double ret = 0.0;
int N = nk * ns;

if (classical) {
if (sum_over_irreducible_points) {

N = kpoint->nk_irred * ns;
int ik_irred;

if (classical) {
#pragma omp parallel for private(ik_irred, ik, is, omega), reduction(+:ret)
for (i = 0; i < N; ++i) {
ik_irred = i / ns;
ik = kpoint->kpoint_irred_all[ik_irred][0].knum;
is = i % ns;

omega = dynamical->eval_phonon[ik][is];

if (omega < 0.0) continue;

ret += Cv_classical(omega, T) * kpoint->weight_k[ik_irred];
}
}
else {
#pragma omp parallel for private(ik_irred, ik, is, omega), reduction(+:ret)
for (i = 0; i < N; ++i) {
ik_irred = i / ns;
ik = kpoint->kpoint_irred_all[ik_irred][0].knum;
is = i % ns;

omega = dynamical->eval_phonon[ik][is];

if (omega < 0.0) continue;

ret += Cv(omega, T) * kpoint->weight_k[ik_irred];
}
}

return ret;
}

if (classical) {
#pragma omp parallel for private(ik, is, omega), reduction(+:ret)
for (i = 0; i < N; ++i) {
ik = i / ns;
is = i % ns;
for (i = 0; i < N; ++i) {
ik = i / ns;
is = i % ns;

omega = dynamical->eval_phonon[ik][is];
if (omega < 0.0) continue;
omega = dynamical->eval_phonon[ik][is];
if (omega < 0.0) continue;

ret += Cv_classical(omega, T);
ret += Cv_classical(omega, T);
}
}
} else {
else {
#pragma omp parallel for private(ik, is, omega), reduction(+:ret)
for (i = 0; i < N; ++i) {
ik = i / ns;
is = i % ns;
for (i = 0; i < N; ++i) {
ik = i / ns;
is = i % ns;

omega = dynamical->eval_phonon[ik][is];
if (omega < 0.0) continue;
omega = dynamical->eval_phonon[ik][is];
if (omega < 0.0) continue;

ret += Cv(omega, T);
ret += Cv(omega, T);
}
}
}

return ret / static_cast<double>(nk);
}
Expand All @@ -110,12 +147,12 @@ double Thermodynamics::Cv_Debye(const double T,
{
unsigned int natmin = system->natmin;
unsigned int i;
double d_theta, theta, theta_max;
double theta, theta_max;
unsigned int ntheta;

double x, y;
double ret;
d_theta = 0.001;
double d_theta = 0.001;

if (TD < eps) {
error->exit("Cv_Debye", "Too small TD");
Expand All @@ -126,28 +163,27 @@ double Thermodynamics::Cv_Debye(const double T,

if (T < eps) {
return 0.0;
} else {
x = TD / T;
theta_max = atan(x);

ret = 0.0;
ntheta = static_cast<unsigned int>(theta_max / d_theta);

for (i = 0; i < ntheta; ++i) {
theta = static_cast<double>(i) * d_theta;
y = tan(theta);
if (y > eps) {
ret += std::pow(y, 4) * std::exp(y)
/ std::pow((std::exp(y) - 1.0) * cos(theta), 2);
}
}
x = TD / T;
theta_max = atan(x);

ret = 0.0;
ntheta = static_cast<unsigned int>(theta_max / d_theta);

for (i = 0; i < ntheta; ++i) {
theta = static_cast<double>(i) * d_theta;
y = tan(theta);
if (y > eps) {
ret += std::pow(y, 4) * std::exp(y)
/ std::pow((std::exp(y) - 1.0) * cos(theta), 2);
}
y = tan(theta_max);
ret += 0.5 * std::pow(y, 4) * std::exp(y)
/ std::pow((std::exp(y) - 1.0) * cos(theta_max), 2);

return 9.0 * static_cast<double>(natmin) * k_Boltzmann
* ret * d_theta / std::pow(x, 3);
}
y = tan(theta_max);
ret += 0.5 * std::pow(y, 4) * std::exp(y)
/ std::pow((std::exp(y) - 1.0) * cos(theta_max), 2);

return 9.0 * static_cast<double>(natmin) * k_Boltzmann
* ret * d_theta / std::pow(x, 3);
}

void Thermodynamics::Debye_T(const double T,
Expand Down Expand Up @@ -181,6 +217,45 @@ double Thermodynamics::internal_energy(const double T)

int N = nk * ns;


if (sum_over_irreducible_points) {


N = kpoint->nk_irred * ns;
int ik_irred;

if (classical) {
#pragma omp parallel for private(ik_irred, ik, is, omega), reduction(+:ret)
for (i = 0; i < N; ++i) {
ik_irred = i / ns;
ik = kpoint->kpoint_irred_all[ik_irred][0].knum;
is = i % ns;
omega = dynamical->eval_phonon[ik][is];

if (omega < eps8) continue;

ret += T_to_Ryd * T * kpoint->weight_k[ik_irred];
}
ret *= 2.0;

} else {
#pragma omp parallel for private(ik_irred, ik, is, omega), reduction(+:ret)
for (i = 0; i < N; ++i) {
ik_irred = i / ns;
ik = kpoint->kpoint_irred_all[ik_irred][0].knum;
is = i % ns;
omega = dynamical->eval_phonon[ik][is];

if (omega < eps8) continue;

ret += omega * coth_T(omega, T) * kpoint->weight_k[ik_irred];
}

}
return ret * 0.5;
}


if (classical) {
#pragma omp parallel for private(ik, is, omega), reduction(+:ret)
for (i = 0; i < N; ++i) {
Expand Down Expand Up @@ -220,6 +295,45 @@ double Thermodynamics::vibrational_entropy(const double T)

int N = nk * ns;


if (sum_over_irreducible_points) {

N = kpoint->nk_irred * ns;
int ik_irred;

if (classical) {
#pragma omp parallel for private(ik_irred, ik, is, omega, x), reduction(+:ret)
for (i = 0; i < N; ++i) {
ik_irred = i / ns;
ik = kpoint->kpoint_irred_all[ik_irred][0].knum;
is = i % ns;
omega = dynamical->eval_phonon[ik][is];

if (omega < eps8 || std::abs(T) < eps) continue;

x = omega / (T * T_to_Ryd);
ret += (std::log(x) - 1.0) * kpoint->weight_k[ik_irred];
}

} else {
#pragma omp parallel for private(ik_irred, ik, is, omega, x), reduction(+:ret)
for (i = 0; i < N; ++i) {
ik_irred = i / ns;
ik = kpoint->kpoint_irred_all[ik_irred][0].knum;
is = i % ns;
omega = dynamical->eval_phonon[ik][is];

if (omega < eps8 || std::abs(T) < eps) continue;

x = omega / (T * T_to_Ryd);
ret += (std::log(1.0 - std::exp(-x)) - x / (std::exp(x) - 1.0)) * kpoint->weight_k[ik_irred];
}
}
return -k_Boltzmann * ret;

}


if (classical) {
#pragma omp parallel for private(ik, is, omega, x), reduction(+:ret)
for (i = 0; i < N; ++i) {
Expand Down Expand Up @@ -261,22 +375,57 @@ double Thermodynamics::free_energy(const double T)

int N = nk * ns;

if (classical) {
#pragma omp parallel for private(ik, is, omega, x), reduction(+:ret)
for (i = 0; i < N; ++i) {
ik = i / ns;
is = i % ns;
omega = dynamical->eval_phonon[ik][is];
if (sum_over_irreducible_points) {

N = kpoint->nk_irred * ns;
int ik_irred;

if (classical) {
#pragma omp parallel for private(ik_irred, ik, is, omega, x), reduction(+:ret)
for (i = 0; i < N; ++i) {
ik_irred = i / ns;
ik = kpoint->kpoint_irred_all[ik_irred][0].knum;
is = i % ns;
omega = dynamical->eval_phonon[ik][is];

if (omega < eps8) continue;

if (std::abs(T) > eps) {
x = omega / (T * T_to_Ryd);
ret += std::log(x) * kpoint->weight_k[ik_irred];
}
}

if (omega < eps8) continue;
return T * T_to_Ryd * ret;

x = omega / (T * T_to_Ryd);
ret += std::log(x);
} else {

#pragma omp parallel for private(ik_irred, ik, is, omega, x), reduction(+:ret)
for (i = 0; i < N; ++i) {
ik_irred = i / ns;
ik = kpoint->kpoint_irred_all[ik_irred][0].knum;
is = i % ns;
omega = dynamical->eval_phonon[ik][is];

if (omega < eps8) continue;

if (std::abs(T) < eps) {
ret += 0.5 * omega * kpoint->weight_k[ik_irred];
}
else {
x = omega / (T * T_to_Ryd);
ret += (0.5 * x + std::log(1.0 - std::exp(-x))) * kpoint->weight_k[ik_irred];
}
}

if (std::abs(T) < eps) return ret;

return T * T_to_Ryd * ret;
}

return T * T_to_Ryd * ret / static_cast<double>(nk);
}

} else {
if (classical) {
#pragma omp parallel for private(ik, is, omega, x), reduction(+:ret)
for (i = 0; i < N; ++i) {
ik = i / ns;
Expand All @@ -285,18 +434,34 @@ double Thermodynamics::free_energy(const double T)

if (omega < eps8) continue;

if (std::abs(T) < eps) {
ret += 0.5 * omega;
} else {
if (std::abs(T) > eps) {
x = omega / (T * T_to_Ryd);
ret += 0.5 * x + std::log(1.0 - std::exp(-x));
ret += std::log(x);
}
}

if (std::abs(T) < eps) return ret / static_cast<double>(nk);

return T * T_to_Ryd * ret / static_cast<double>(nk);

}
#pragma omp parallel for private(ik, is, omega, x), reduction(+:ret)
for (i = 0; i < N; ++i) {
ik = i / ns;
is = i % ns;
omega = dynamical->eval_phonon[ik][is];

if (omega < eps8) continue;

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

if (std::abs(T) < eps) return ret / static_cast<double>(nk);

return T * T_to_Ryd * ret / static_cast<double>(nk);
}

double Thermodynamics::disp2_avg(const double T,
Expand Down
1 change: 1 addition & 0 deletions anphon/thermodynamics.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ namespace PHON_NS

double T_to_Ryd;
bool classical;
bool sum_over_irreducible_points;
void setup();

double Cv(double, double);
Expand Down

0 comments on commit 45a6c63

Please sign in to comment.