Skip to content

Commit

Permalink
Implement default values of ISOFACT for most of the elements
Browse files Browse the repository at this point in the history
  • Loading branch information
@ttadano
ttadano committed Feb 23, 2018
1 parent 341fc62 commit 01ddb34
Show file tree
Hide file tree
Showing 4 changed files with 130 additions and 98 deletions.
31 changes: 30 additions & 1 deletion anphon/isotope.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#include "isotope.h"
#include "constants.h"
#include "dynamical.h"
#include "error.h"
#include "integration.h"
#include "kpoint.h"
#include "memory.h"
Expand All @@ -29,7 +30,6 @@ Isotope::Isotope(PHON *phon): Pointers(phon)
Isotope::~Isotope()
{
deallocate_variables();

};

void Isotope::set_default_variables()
Expand Down Expand Up @@ -59,6 +59,13 @@ void Isotope::setup_isotope_scattering()

if (include_isotope) {

if (mympi->my_rank == 0) {
if (!isotope_factor) {
memory->allocate(isotope_factor, nkd);
set_isotope_factor_from_database(nkd, system->symbol_kd, isotope_factor);
}
}

if (mympi->my_rank > 0) {
memory->allocate(isotope_factor, nkd);
}
Expand Down Expand Up @@ -238,3 +245,25 @@ void Isotope::calc_isotope_selfenergy_all()
}
}
}

void Isotope::set_isotope_factor_from_database(const int nkd,
const std::string *symbol_in,
double *isofact_out)
{
double isofact_tmp;
int atom_number;

for (int i = 0; i < nkd; ++i) {
atom_number = system->get_atomic_number_by_name(symbol_in[i]);
if (atom_number >= isotope_factors.size() || (atom_number == -1)) {
error->exit("set_isotope_factor_from_database",
"The isotope factor for the given element doesn't exist in the database.\nTherefore, please input ISOFACT manually.");
}
isofact_tmp = isotope_factors[atom_number];
if (isofact_tmp < -0.5) {
error->exit("set_isotope_factor_from_database",
"One of the elements in the KD-tag is unstable. \nTherefore, please input ISOFACT manually.");
}
isofact_out[i] = isofact_tmp;
}
}
24 changes: 23 additions & 1 deletion anphon/isotope.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
*/

#include "pointers.h"

#include <vector>

namespace PHON_NS
{
Expand Down Expand Up @@ -39,5 +39,27 @@ namespace PHON_NS
int,
double,
double &);
void set_isotope_factor_from_database(const int,
const std::string *,
double *);

std::vector<double> isotope_factors{
0.00011460742533846188, 8.141017510851699e-08, 0.001458820047621947, 0.0, 0.0013539153245004177,
7.387230454537788e-05, 1.8376657914224594e-05, 3.3588025065260784e-05, 0.0, 0.0008278312160539562, 0.0,
0.0007398827105109379, 0.0, 0.00020070043953754622, 0.0, 0.00016512811739588583, 0.0005827012790438193,
3.4803742597617525e-05, 0.0001640001019257936, 0.00029756377957496945, 0.0, 0.00028645556480110174,
9.548320927761995e-07, 0.0001328704171052577, 0.0, 8.244411853712163e-05, 0.0, 0.00043070442089852476,
0.00021093312802220393, 0.0005931533235650079, 0.00019712731536263687, 0.000586966208421238, 0.0,
0.0004627901461335774, 0.00015627716697599108, 0.00024880742737704476, 0.00010969638987469167,
6.099700015374489e-05, 0.0, 0.00034262903024295327, 0.0, 0.0005961083777486782, -1, 0.0004066661504740231,
0.0, 0.0003094784411091063, 8.579847704787673e-05, 0.0002716036180261363, 1.245588189674909e-05,
0.0003340852797872776, 6.607553852631361e-05, 0.0002839333030058624, 0.0, 0.0002675566535085368, 0.0,
6.237013178021676e-05, 4.5917491111023726e-08, 2.2495590932891925e-05, 0.0, 0.00023237187990100274, -1,
0.000334685954430736, 4.3279441126609935e-05, 0.0001276749037273739, 0.0, 5.198070714285335e-05, 0.0,
7.23248017569606e-05, 0.0, 8.556028002982757e-05, 8.300794202558322e-07, 5.253850496170609e-05,
3.6715121208243084e-09, 6.966807117351981e-05, 2.7084982818795603e-05, 7.452354225251159e-05,
2.5378700157091918e-05, 3.4285145177491544e-05, 0.0, 6.525193204276608e-05, 1.9964351041965618e-05,
1.9437780365209887e-05, 0.0, -1, -1, -1, -1, -1, -1, 0.0, 0.0, 1.1564592331193284e-06, -1, -1
}; // Precalculated isotope factors. For unstable elements, the value is set to -1.
};
}
152 changes: 66 additions & 86 deletions anphon/parsephon.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -327,19 +327,12 @@ void Input::parse_scph_vars()
// Read input parameters in the &scph-field.

int i;
unsigned int maxiter;
unsigned int ialgo_scph;
double tolerance_scph;
double mixalpha;
bool restart_scph;
bool selenergy_offdiagonal;
bool lower_temp, warm_start;

struct stat st;
std::string file_dymat;
std::string str_tmp;
std::vector<std::string> input_list{
"KMESH_SCPH", "KMESH_INTERPOLATE", "MIXALPHA", "MAXITER", "RESTART_SCPH", "IALGO", "SELF_OFFDIAG", "TOL_SCPH",
"KMESH_SCPH", "KMESH_INTERPOLATE", "MIXALPHA", "MAXITER",
"RESTART_SCPH", "IALGO", "SELF_OFFDIAG", "TOL_SCPH",
"LOWER_TEMP", "WARMSTART"
};
std::vector<std::string> no_defaults{"KMESH_SCPH", "KMESH_INTERPOLATE"};
Expand All @@ -353,40 +346,31 @@ void Input::parse_scph_vars()
}

get_var_dict(input_list, scph_var_dict);
//#if _USE_BOOST
// boost::split(no_defaults, str_no_defaults, boost::is_space());
//#else
// no_defaults = my_split(str_no_defaults, ' ');
//#endif

for (auto it = no_defaults.begin(); it != no_defaults.end(); ++it) {
if (scph_var_dict.find(*it) == scph_var_dict.end()) {

for (auto &no_default : no_defaults) {
if (scph_var_dict.find(no_default) == scph_var_dict.end()) {
error->exit("parse_general_vars",
"The following variable is not found in &scph input region: ",
(*it).c_str());
no_default.c_str());
}
}

file_dymat = this->job_title + ".scph_dymat";
std::string file_dymat = this->job_title + ".scph_dymat";

// Default values

tolerance_scph = 1.0e-10;
maxiter = 1000;
mixalpha = 0.1;
selenergy_offdiagonal = true;
ialgo_scph = 0;
lower_temp = true;
warm_start = true;
double tolerance_scph = 1.0e-10;
unsigned int maxiter = 1000;
double mixalpha = 0.1;
bool selenergy_offdiagonal = true;
unsigned int ialgo_scph = 0;
bool lower_temp = true;
bool warm_start = true;

// if file_dymat exists in the current directory,
// restart mode will be automatically turned on for SCPH calculations.

if (stat(file_dymat.c_str(), &st) == 0) {
restart_scph = true;
} else {
restart_scph = false;
}
bool restart_scph = stat(file_dymat.c_str(), &st) == 0;

// Assign given values

Expand Down Expand Up @@ -474,66 +458,53 @@ void Input::parse_analysis_vars(const bool use_default_values)
int i;

std::vector<std::string> input_list{
"PRINTEVEC", "PRINTXSF", "PRINTVEL", "QUARTIC", "KS_INPUT", "ATOMPROJ", "REALPART", "ISOTOPE", "ISOFACT",
"FSTATE_W", "FSTATE_K", "PRIMTMSD", "PDOS", "TDOS", "GRUNEISEN", "NEWFCS", "DELTA_A", "ANIME", "ANIME_CELLSIZE",
"ANIME_FORMAT", "SPS", "PRINTV3", "PRINTPR", "FC2_EWALD", "KAPPA_SPEC", "SELF_W"
"PRINTEVEC", "PRINTXSF", "PRINTVEL", "QUARTIC", "KS_INPUT",
"ATOMPROJ", "REALPART", "ISOTOPE", "ISOFACT",
"FSTATE_W", "FSTATE_K", "PRIMTMSD", "PDOS", "TDOS",
"GRUNEISEN", "NEWFCS", "DELTA_A", "ANIME", "ANIME_CELLSIZE",
"ANIME_FORMAT", "SPS", "PRINTV3", "PRINTPR", "FC2_EWALD",
"KAPPA_SPEC", "SELF_W"
};

bool fstate_omega, fstate_k;
bool ks_analyze_mode, atom_project_mode, calc_realpart;
bool print_vel, print_evec, print_msd;
bool projected_dos, print_gruneisen, print_newfcs;
bool two_phonon_dos;
bool print_xsf, print_anime;
bool print_V3, participation_ratio;
bool print_fc2_ewald;
bool print_self_consistent_fc2;
bool bubble_omega;

int quartic_mode;
int include_isotope;
int scattering_phase_space;
int calculate_kappa_spec;
unsigned int cellsize[3];

double delta_a;
double *isotope_factor;
double *isotope_factor = nullptr;
std::string ks_input, anime_format;
std::map<std::string, std::string> analysis_var_dict;
std::vector<std::string> isofact_v, anime_kpoint, anime_cellsize;

// Default values

print_xsf = false;
print_anime = false;
bool print_xsf = false;
bool print_anime = false;

print_vel = false;
print_evec = false;
print_msd = false;
bool print_vel = false;
bool print_evec = false;
bool print_msd = false;

projected_dos = false;
two_phonon_dos = false;
scattering_phase_space = 0;
print_gruneisen = false;
print_newfcs = false;
print_V3 = false;
participation_ratio = false;
bool projected_dos = false;
bool two_phonon_dos = false;
int scattering_phase_space = 0;
bool print_gruneisen = false;
bool print_newfcs = false;
bool print_V3 = false;
bool participation_ratio = false;

delta_a = 0.001;
double delta_a = 0.001;

quartic_mode = 0;
ks_analyze_mode = false;
atom_project_mode = false;
calc_realpart = false;
include_isotope = 0;
fstate_omega = false;
fstate_k = false;
bubble_omega = false;
int quartic_mode = 0;
bool ks_analyze_mode = false;
bool atom_project_mode = false;
bool calc_realpart = false;
int include_isotope = 0;
bool fstate_omega = false;
bool fstate_k = false;
bool bubble_omega = false;

calculate_kappa_spec = 0;
int calculate_kappa_spec = 0;

print_fc2_ewald = false;
print_self_consistent_fc2 = false;
bool print_fc2_ewald = false;
bool print_self_consistent_fc2 = false;


// Assign values to variables
Expand Down Expand Up @@ -575,17 +546,22 @@ void Input::parse_analysis_vars(const bool use_default_values)
}

if (include_isotope) {
split_str_by_space(analysis_var_dict["ISOFACT"], isofact_v);

if (isofact_v.size() != system->nkd) {
error->exit("parse_analysis_vars",
"The number of entries for ISOFACT is inconsistent with NKD");
} else {
memory->allocate(isotope_factor, system->nkd);
for (i = 0; i < system->nkd; ++i) {
isotope_factor[i] = my_cast<double>(isofact_v[i]);
if (!analysis_var_dict["ISOFACT"].empty()) {
split_str_by_space(analysis_var_dict["ISOFACT"], isofact_v);

if (isofact_v.size() != system->nkd) {
error->exit("parse_analysis_vars",
"The number of entries for ISOFACT is inconsistent with NKD");
}
else {
memory->allocate(isotope_factor, system->nkd);
for (i = 0; i < system->nkd; ++i) {
isotope_factor[i] = my_cast<double>(isofact_v[i]);
}
}
}

}

if (print_anime) {
Expand Down Expand Up @@ -669,10 +645,14 @@ void Input::parse_analysis_vars(const bool use_default_values)
ewald->print_fc2_ewald = print_fc2_ewald;

if (include_isotope) {
memory->allocate(isotope->isotope_factor, system->nkd);
for (i = 0; i < system->nkd; ++i) {
isotope->isotope_factor[i] = isotope_factor[i];
if (!analysis_var_dict["ISOFACT"].empty()) {
memory->allocate(isotope->isotope_factor, system->nkd);
for (i = 0; i < system->nkd; ++i) {
isotope->isotope_factor[i] = isotope_factor[i];
}
}
}
if (isotope_factor) {
memory->deallocate(isotope_factor);
}

Expand Down
21 changes: 11 additions & 10 deletions anphon/write_phonons.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,9 +43,7 @@ Writes::Writes(PHON *phon): Pointers(phon)
Ry_to_kayser = Hz_to_kayser / time_ry;
};

Writes::~Writes()
{
};
Writes::~Writes() {};

void Writes::write_input_vars()
{
Expand All @@ -69,8 +67,10 @@ void Writes::write_input_vars()
}
std::cout << std::endl;
std::cout << " MASS = ";
for (i = 0; i < system->nkd; ++i) {
std::cout << std::setw(10) << system->mass_kd[i];
if (system->mass_kd) {
for (i = 0; i < system->nkd; ++i) {
std::cout << std::setw(10) << system->mass_kd[i];
}
}
std::cout << std::endl;
std::cout << " NSYM = " << symmetry->nsym << "; TOLERANCE = " << symmetry->tolerance;
Expand Down Expand Up @@ -170,9 +170,11 @@ void Writes::write_input_vars()
std::cout << " ISOTOPE = " << isotope->include_isotope << std::endl;
if (isotope->include_isotope) {
std::cout << " ISOFACT = ";
for (i = 0; i < system->nkd; ++i) {
std::cout << std::scientific
<< std::setw(13) << isotope->isotope_factor[i];
if (isotope->isotope_factor) {
for (i = 0; i < system->nkd; ++i) {
std::cout << std::scientific
<< std::setw(13) << isotope->isotope_factor[i];
}
}
std::cout << std::endl;
}
Expand All @@ -187,8 +189,7 @@ void Writes::write_input_vars()
// std::cout << " FSTATE_K = " << relaxation->calc_fstate_k << std::endl;

} else if (phon->mode == "SCPH") {


// Do nothing
} else {
error->exit("write_input_vars", "This cannot happen");
}
Expand Down

0 comments on commit 01ddb34

Please sign in to comment.