Skip to content

Commit

Permalink
Restore the CLASSICAL option
Browse files Browse the repository at this point in the history
  • Loading branch information
@ttadano
ttadano committed Nov 20, 2017
1 parent 567923f commit eedde28
Show file tree
Hide file tree
Showing 12 changed files with 323 additions and 94 deletions.
26 changes: 21 additions & 5 deletions anphon/conductivity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -351,6 +351,7 @@ void Conductivity::compute_kappa()
std::string file_kl;
std::ofstream ofs_kl;
double vv_tmp;
double damp_tmp;
int ieq;

double **lifetime;
Expand All @@ -370,8 +371,12 @@ void Conductivity::compute_kappa()
}
} else {
for (i = 0; i < ntemp; ++i) {
lifetime[iks][i] = 1.0e+12 * time_ry * 0.5
/ (damping3[iks][i] + isotope->gamma_isotope[iks / ns][snum]);
damp_tmp = damping3[iks][i] + isotope->gamma_isotope[iks / ns][snum];
if (damp_tmp > 1.0e-100) {
lifetime[iks][i] = 1.0e+12 * time_ry * 0.5 / damp_tmp;
} else {
lifetime[iks][i] = 0.0;
}
}
}
}
Expand All @@ -384,7 +389,12 @@ void Conductivity::compute_kappa()
}
} else {
for (i = 0; i < ntemp; ++i) {
lifetime[iks][i] = 1.0e+12 * time_ry * 0.5 / damping3[iks][i];
damp_tmp = damping3[iks][i];
if (damp_tmp > 1.0e-100) {
lifetime[iks][i] = 1.0e+12 * time_ry * 0.5 / damp_tmp;
} else {
lifetime[iks][i] = 0.0;
}
}
}
}
Expand Down Expand Up @@ -417,8 +427,14 @@ void Conductivity::compute_kappa()
vv_tmp += vel[ktmp][is][j] * vel[ktmp][is][k];
}

kappa_mode[i][3 * j + k][is][ik] = thermodynamics->Cv(omega, Temperature[i])
* vv_tmp * lifetime[ns * ik + is][i];
if (thermodynamics->classical) {
kappa_mode[i][3 * j + k][is][ik] = thermodynamics->Cv_classical(omega, Temperature[i])
* vv_tmp * lifetime[ns * ik + is][i];
} else {
kappa_mode[i][3 * j + k][is][ik] = thermodynamics->Cv(omega, Temperature[i])
* vv_tmp * lifetime[ns * ik + is][i];
}


// Convert to SI unit
kappa_mode[i][3 * j + k][is][ik] *= factor_toSI;
Expand Down
8 changes: 6 additions & 2 deletions anphon/parsephon.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
#include "integration.h"
#include "scph.h"
#include "ewald.h"
#include "thermodynamics.h"
#include <boost/lexical_cast.hpp>

#include <boost/algorithm/string.hpp>
Expand Down Expand Up @@ -118,6 +119,7 @@ void Input::parse_general_vars()
bool sym_time_reversal, use_triplet_symmetry;
bool selenergy_offdiagonal;
bool update_fc2;
bool classical;

struct stat st;
std::string prefix, mode, fcsinfo, fc2info;
Expand All @@ -126,7 +128,7 @@ void Input::parse_general_vars()
std::string str_tmp;
std::string str_allowed_list = "PREFIX MODE NSYM TOLERANCE PRINTSYM FCSXML FC2XML TMIN TMAX DT \
NBANDS NONANALYTIC BORNINFO NA_SIGMA ISMEAR EPSILON EMIN EMAX DELTA_E \
RESTART TREVSYM NKD KD MASS TRISYM PREC_EWALD";
RESTART TREVSYM NKD KD MASS TRISYM PREC_EWALD CLASSICAL";
std::string str_no_defaults = "PREFIX MODE FCSXML NKD KD MASS";
std::vector<std::string> no_defaults;
std::vector<std::string> kdname_v, masskd_v;
Expand Down Expand Up @@ -204,6 +206,7 @@ void Input::parse_general_vars()
printsymmetry = false;
sym_time_reversal = false;
use_triplet_symmetry = true;
classical = false;

prec_ewald = 1.0e-12;

Expand Down Expand Up @@ -249,6 +252,7 @@ void Input::parse_general_vars()
assign_val(ismear, "ISMEAR", general_var_dict);
assign_val(epsilon, "EPSILON", general_var_dict);
assign_val(na_sigma, "NA_SIGMA", general_var_dict);
assign_val(classical, "CLASSICAL", general_var_dict);

assign_val(use_triplet_symmetry, "TRISYM", general_var_dict);

Expand Down Expand Up @@ -319,7 +323,7 @@ void Input::parse_general_vars()
}
fcs_phonon->file_fc2 = fc2info;
fcs_phonon->update_fc2 = update_fc2;

thermodynamics->classical = classical;
integration->ismear = ismear;
relaxation->use_triplet_symmetry = use_triplet_symmetry;

Expand Down
16 changes: 12 additions & 4 deletions anphon/phonon_dos.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -813,10 +813,18 @@ void Dos::calc_scattering_phase_space_with_Bose_mode(const unsigned int nk,

if (omega1 < eps12 || omega2 < eps12) continue;

f1 = thermodynamics->fB(omega1, temp);
f2 = thermodynamics->fB(omega2, temp);
n1 = f1 + f2 + 1.0;
n2 = f1 - f2;
if (thermodynamics->classical) {
f1 = thermodynamics->fC(omega1, temp);
f2 = thermodynamics->fC(omega2, temp);
n1 = f1 + f2;
n2 = f1 - f2;
} else {
f1 = thermodynamics->fB(omega1, temp);
f2 = thermodynamics->fB(omega2, temp);
n1 = f1 + f2 + 1.0;
n2 = f1 - f2;
}


ret1 += delta_arr[k1][ib][0] * n1;
ret2 += -delta_arr[k1][ib][1] * n2;
Expand Down
5 changes: 5 additions & 0 deletions anphon/phonons.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,9 +158,14 @@ void PHON::setup_base()
fcs_phonon->setup(mode);
dynamical->setup_dynamical(mode);
dos->setup();
thermodynamics->setup();
ewald->init();
if (mympi->my_rank == 0) {
std::cout << " Now, move on to phonon calculations." << std::endl;
if (thermodynamics->classical) {
std::cout << std::endl;
std::cout << " CLASSICAL = 1: Classical approximations will be used for all thermodynamical quantities." << std::endl;
}
}
}

Expand Down
99 changes: 74 additions & 25 deletions anphon/relaxation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -921,15 +921,24 @@ void Relaxation::calc_damping_smearing(const unsigned int N,
for (is = 0; is < ns; ++is) {

omega_inner[0] = dynamical->eval_phonon[k1][is];
f1 = thermodynamics->fB(omega_inner[0], T_tmp);

for (js = 0; js < ns; ++js) {

omega_inner[1] = dynamical->eval_phonon[k2][js];
f2 = thermodynamics->fB(omega_inner[1], T_tmp);

n1 = f1 + f2 + 1.0;
n2 = f1 - f2;
if (thermodynamics->classical) {
f1 = thermodynamics->fC(omega_inner[0], T_tmp);
f2 = thermodynamics->fC(omega_inner[1], T_tmp);

n1 = f1 + f2;
n2 = f1 - f2;
} else {
f1 = thermodynamics->fB(omega_inner[0], T_tmp);
f2 = thermodynamics->fB(omega_inner[1], T_tmp);

n1 = f1 + f2 + 1.0;
n2 = f1 - f2;
}

ret_tmp += v3_arr[ik][ns * is + js]
* (n1 * delta_arr[ik][ns * is + js][0]
Expand Down Expand Up @@ -1086,15 +1095,24 @@ void Relaxation::calc_damping_tetrahedron(const unsigned int N,
for (is = 0; is < ns; ++is) {

omega_inner[0] = dynamical->eval_phonon[k1][is];
f1 = thermodynamics->fB(omega_inner[0], T_tmp);

for (js = 0; js < ns; ++js) {

omega_inner[1] = dynamical->eval_phonon[k2][js];
f2 = thermodynamics->fB(omega_inner[1], T_tmp);

n1 = f1 + f2 + 1.0;
n2 = f1 - f2;
if (thermodynamics->classical) {
f1 = thermodynamics->fC(omega_inner[0], T_tmp);
f2 = thermodynamics->fC(omega_inner[1], T_tmp);

n1 = f1 + f2;
n2 = f1 - f1;
} else {
f1 = thermodynamics->fB(omega_inner[0], T_tmp);
f2 = thermodynamics->fB(omega_inner[1], T_tmp);

n1 = f1 + f2 + 1.0;
n2 = f1 - f2;
}

ret_tmp += v3_arr[ik][ns * is + js]
* (n1 * delta_arr[ik][ns * is + js][0]
Expand Down Expand Up @@ -1179,11 +1197,18 @@ void Relaxation::calc_frequency_resolved_final_state(const unsigned int N,
for (i = 0; i < N; ++i) {
T_tmp = T[i];

f1 = thermodynamics->fB(omega_inner[0], T_tmp);
f2 = thermodynamics->fB(omega_inner[1], T_tmp);
n1 = f1 + f2 + 1.0;
n2 = f1 - f2;

if (thermodynamics->classical) {
f1 = thermodynamics->fC(omega_inner[0], T_tmp);
f2 = thermodynamics->fC(omega_inner[1], T_tmp);
n1 = f1 + f2;
n2 = f1 - f2;
} else {
f1 = thermodynamics->fB(omega_inner[0], T_tmp);
f2 = thermodynamics->fB(omega_inner[1], T_tmp);
n1 = f1 + f2 + 1.0;
n2 = f1 - f2;
}

if (integration->ismear == 0) {
prod_tmp[0] = n1
* (delta_lorentz(omega0 - omega_inner[0] - omega_inner[1], epsilon)
Expand Down Expand Up @@ -1388,11 +1413,19 @@ void Relaxation::calc_frequency_resolved_final_state_tetrahedron(const unsigned
#endif
for (i = 0; i < N; ++i) {

f1 = thermodynamics->fB(omega_inner[0], T[i]);
f2 = thermodynamics->fB(omega_inner[1], T[i]);
if (thermodynamics->classical) {
f1 = thermodynamics->fC(omega_inner[0], T[i]);
f2 = thermodynamics->fC(omega_inner[1], T[i]);

n1 = f1 + f2 + 1.0;
n2 = f1 - f2;
n1 = f1 + f2;
n2 = f1 - f2;
} else {
f1 = thermodynamics->fB(omega_inner[0], T[i]);
f2 = thermodynamics->fB(omega_inner[1], T[i]);

n1 = f1 + f2 + 1.0;
n2 = f1 - f2;
}

prod_tmp[0] = v3_tmp * n1 * delta_arr[ik][ns * is + js][0];
prod_tmp[1] = -v3_tmp * n2 * delta_arr[ik][ns * is + js][1];
Expand Down Expand Up @@ -2498,10 +2531,18 @@ void Relaxation::print_momentum_resolved_final_state(const unsigned int NT,
for (iT = 0; iT < NT; ++iT) {
T_tmp = T_arr[iT];

f1 = thermodynamics->fB(eval[1][is], T_tmp);
f2 = thermodynamics->fB(eval[2][js], T_tmp);
n1 = f1 + f2 + 1.0;
n2 = f1 - f2;
if (thermodynamics->classical) {
f1 = thermodynamics->fC(eval[1][is], T_tmp);
f2 = thermodynamics->fC(eval[2][js], T_tmp);
n1 = f1 + f2;
n2 = f1 - f2;
} else {
f1 = thermodynamics->fB(eval[1][is], T_tmp);
f2 = thermodynamics->fB(eval[2][js], T_tmp);
n1 = f1 + f2 + 1.0;
n2 = f1 - f2;
}


if (selection_type == 0) {
gamma_k[k][iT] += V3norm * n1;
Expand Down Expand Up @@ -3068,10 +3109,18 @@ void Relaxation::calc_self3omega_tetrahedron(const double Temp,

omega_inner[0] = eval[k1][is];
omega_inner[1] = eval[k2][js];
f1 = thermodynamics->fB(omega_inner[0], Temp);
f2 = thermodynamics->fB(omega_inner[1], Temp);
n1 = f1 + f2 + 1.0;
n2 = f1 - f2;
if (thermodynamics->classical) {
f1 = thermodynamics->fC(omega_inner[0], Temp);
f2 = thermodynamics->fC(omega_inner[1], Temp);
n1 = f1 + f2;
n2 = f1 - f2;
} else {
f1 = thermodynamics->fB(omega_inner[0], Temp);
f2 = thermodynamics->fB(omega_inner[1], Temp);
n1 = f1 + f2 + 1.0;
n2 = f1 - f2;
}

//#pragma omp critical
ret_private[nomega * ithread + iomega]
+= v3_arr[ik][ib] * (n1 * weight_tetra[0][ik] - 2.0 * n2 * weight_tetra[1][ik]);
Expand Down
23 changes: 17 additions & 6 deletions anphon/scph.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@ void Scph::exec_scph()
write_scph_bands(eval_anharm);
} else if (kpoint->kpoint_mode == 2) {
write_scph_dos(eval_anharm);
write_scph_thermodynamics(eval_anharm);
// write_scph_thermodynamics(eval_anharm);
if (writes->print_msd) {
write_scph_msd(eval_anharm, evec_anharm);
}
Expand Down Expand Up @@ -2529,8 +2529,14 @@ void Scph::compute_anharmonic_frequency(std::complex<double> ***v4_array_all,
if (std::abs(omega1) < eps8) {
Qmat(is, is) = complex_zero;
} else {
n1 = thermodynamics->fB(omega1, T_in);
Qmat(is, is) = std::complex<double>((2.0 * n1 + 1.0) / omega1, 0.0);
// Note that the missing factor 2 in the denominator of Qmat is
// already considered in the v4_array_all.
if (thermodynamics->classical) {
Qmat(is, is) = std::complex<double>(2.0 * T_in * thermodynamics->T_to_Ryd / (omega1 * omega1), 0.0);
} else {
n1 = thermodynamics->fB(omega1, T_in);
Qmat(is, is) = std::complex<double>((2.0 * n1 + 1.0) / omega1, 0.0);
}
}
}

Expand Down Expand Up @@ -3067,6 +3073,7 @@ void Scph::write_scph_dos(double ***eval)

void Scph::write_scph_thermodynamics(double ***eval)
{
// This function is incorrect. Don't use it.
int i;
unsigned int iT;
unsigned int ik, is;
Expand Down Expand Up @@ -3171,9 +3178,13 @@ void Scph::write_scph_msd(double ***eval,
omega = eval[iT][ik][is];

if (omega < eps8) continue;

tmp += std::norm(evec[iT][ik][is][i])
* (thermodynamics->fB(omega, temp) + 0.5) / omega;
if (thermodynamics->classical) {
tmp += std::norm(evec[iT][ik][is][i])
* thermodynamics->fC(omega, temp) / omega;
} else {
tmp += std::norm(evec[iT][ik][is][i])
* (thermodynamics->fB(omega, temp) + 0.5) / omega;
}
}
}
msd[iT][i] = tmp * std::pow(Bohr_in_Angstrom, 2.0)
Expand Down

0 comments on commit eedde28

Please sign in to comment.