Skip to content

Commit

Permalink
Testing the connection of band dispersion using the similarity of eig…
Browse files Browse the repository at this point in the history 
…envectors
  • Loading branch information
@ttadano
ttadano committed Nov 28, 2017
1 parent 4d3454d commit 3f41331
Show file tree
Hide file tree
Showing 4 changed files with 162 additions and 12 deletions.
113 changes: 108 additions & 5 deletions anphon/dynamical.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,12 +28,14 @@
#include "phonon_dos.h"
#include "gruneisen.h"
#include "ewald.h"
#include <numeric>


using namespace PHON_NS;

Dynamical::Dynamical(PHON *phon): Pointers(phon)
{
index_bconnect = nullptr;
}

Dynamical::~Dynamical()
Expand Down Expand Up @@ -114,6 +116,7 @@ void Dynamical::setup_dynamical(std::string mode)

MPI_Bcast(&eigenvectors, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD);
MPI_Bcast(&nonanalytic, 1, MPI_UNSIGNED, 0, MPI_COMM_WORLD);
MPI_Bcast(&band_connection, 1, MPI_UNSIGNED, 0, MPI_COMM_WORLD);

if (nonanalytic) {
memory->allocate(borncharge, system->natmin, 3, 3);
Expand Down Expand Up @@ -148,6 +151,10 @@ void PHON_NS::Dynamical::finish_dynamical()
if (nonanalytic) {
memory->deallocate(borncharge);
}

if (index_bconnect) {
memory->deallocate(index_bconnect);
}
}


Expand Down Expand Up @@ -741,6 +748,11 @@ void Dynamical::diagonalize_dynamical_all()
}
}

if (band_connection > 0 && kpoint->kpoint_mode == 1) {
memory->allocate(index_bconnect, nk, neval);
connect_band_by_eigen_similarity(evec_phonon, index_bconnect);
}

if (mympi->my_rank == 0) {
std::cout << "done!" << std::endl;
}
Expand Down Expand Up @@ -852,8 +864,8 @@ void Dynamical::load_born()

for (j = 0; j < 3; ++j) {
for (k = 0; k < 3; ++k) {
std::cout << std::setw(15) << std::fixed
<< std::setprecision(6) << borncharge[i][j][k];
std::cout << std::setw(15) << std::fixed
<< std::setprecision(6) << borncharge[i][j][k];
}
std::cout << std::endl;
}
Expand Down Expand Up @@ -930,7 +942,7 @@ void Dynamical::load_born()
}
}
}

for (iat = 0; iat < system->natmin; ++iat) {
for (i = 0; i < 3; ++i) {
for (j = 0; j < 3; ++j) {
Expand Down Expand Up @@ -963,14 +975,14 @@ void Dynamical::load_born()
}
}
memory->deallocate(born_sym);

if (diff_sym > eps8 || res > eps10) {
std::cout << std::endl;
std::cout << " Symmetrized Born effective charge tensor in Cartesian coordinate." << std::endl;
for (i = 0; i < system->natmin; ++i) {
std::cout << " Atom" << std::setw(5) << i + 1 << "("
<< std::setw(3) << system->symbol_kd[system->kd[system->map_p2s[i][0]]] << ") :" << std::endl;

for (j = 0; j < 3; ++j) {
for (k = 0; k < 3; ++k) {
std::cout << std::setw(15) << borncharge[i][j][k];
Expand Down Expand Up @@ -1055,3 +1067,94 @@ void Dynamical::calc_atomic_participation_ratio(std::complex<double> *evec, doub
for (iat = 0; iat < natmin; ++iat)
ret[iat] /= std::sqrt(static_cast<double>(natmin) * sum);
}


void Dynamical::connect_band_by_eigen_similarity(std::complex<double> ***evec,
int **index_sorted)
{
int ik, is, js, ks;
unsigned int nk = kpoint->nk;
unsigned int ns = neval;
int loc;
std::vector<int> index;
std::complex<double> **evec_tmp;
std::vector<std::vector<double>> abs_similarity;
std::complex<double> dprod;
std::vector<int> found;

memory->allocate(evec_tmp, ns, ns);

for (ik = 0; ik < nk; ++ik) {
for (is = 0; is < ns; ++is) {
index_sorted[ik][is] = 0;
}
}

index.resize(ns);
found.resize(ns);
abs_similarity.resize(ns);
for (is = 0; is < ns; ++is) {
abs_similarity[is].resize(ns);
}

for (int i = 0; i < ns; ++i) index[i] = i;

for (ik = 0; ik < nk; ++ik) {

if (ik == 0) {
for (is = 0; is < ns; ++is) {
for (js = 0; js < ns; ++js) {
if (is == js) {
abs_similarity[is][js] = 1.0;
} else {
abs_similarity[is][js] = 0.0;
}
}
}
} else {
//#ifdef _OPENMP
//#pragma omp parallel for private(js, ks, dprod)
//#endif
for (is = 0; is < ns; ++is) {
for (js = 0; js < ns; ++js) {
dprod = std::complex<double>(0.0, 0.0);
for (ks = 0; ks < ns; ++ks) {
dprod += std::conj(evec[ik][is][ks]) * evec_tmp[js][ks];
}
abs_similarity[is][js] = std::abs(dprod);
}
}
}

for (auto &v : found) v = 0;

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

for (auto it = abs_similarity[is].begin(); it != abs_similarity[is].end(); ++it) {
std::cout << (*it) << std::endl;
}
std::cout << std::endl;
// Argsort abs_similarity[is] (use C++11 lambda)
iota(index.begin(), index.end(), 0);
std::sort(index.begin(), index.end(),
[&abs_similarity, is](int i1, int i2) {
return abs_similarity[is][i1] > abs_similarity[is][i2];
});

loc = index[0];
index_sorted[ik][loc] = is;
found[loc] = 1;
for (js = 0; js < ns; ++js) abs_similarity[js][loc] = -1.0;
for (js = 0; js < ns; ++js) {
evec_tmp[loc][js] = evec[ik][is][js];
}
}

if (std::any_of(found.begin(), found.end(), [](int i1) { return i1 == 0; })) {
error->exit("connect_band_by_eigen_similarity",
"Could not identify the connection.");
}

}
memory->deallocate(evec_tmp);
}
3 changes: 3 additions & 0 deletions anphon/dynamical.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,11 +53,13 @@ namespace PHON_NS
bool print_eigenvectors;
unsigned int nonanalytic;
bool participation_ratio;
unsigned int band_connection;

std::string file_born;
double na_sigma;

double **eval_phonon;
int **index_bconnect;
std::complex<double> ***evec_phonon;
double dielec[3][3];
double ***borncharge;
Expand Down Expand Up @@ -101,6 +103,7 @@ namespace PHON_NS
void prepare_mindist_list(std::vector<int> **);
void calc_atomic_participation_ratio(std::complex<double> *, double *);
double distance(double *, double *);
void connect_band_by_eigen_similarity(std::complex<double> ***, int **);

// void calc_analytic_k(double *, double ****, std::complex<double> **);
// void modify_eigenvectors_sym();
Expand Down
11 changes: 9 additions & 2 deletions anphon/parsephon.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,6 +120,7 @@ void Input::parse_general_vars()
bool selenergy_offdiagonal;
bool update_fc2;
bool classical;
unsigned int band_connection;

struct stat st;
std::string prefix, mode, fcsinfo, fc2info;
Expand All @@ -128,7 +129,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 CLASSICAL";
RESTART TREVSYM NKD KD MASS TRISYM PREC_EWALD CLASSICAL BCONNECT";
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 @@ -207,6 +208,7 @@ void Input::parse_general_vars()
sym_time_reversal = false;
use_triplet_symmetry = true;
classical = false;
band_connection = 0;

prec_ewald = 1.0e-12;

Expand Down Expand Up @@ -253,9 +255,13 @@ void Input::parse_general_vars()
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(band_connection, "BCONNECT", general_var_dict);
assign_val(use_triplet_symmetry, "TRISYM", general_var_dict);

if (band_connection < 0 || band_connection > 2) {
error->exit("parse_general_vars", "BCONNECT-tag can take 0, 1, or 2.");
}

if (nonanalytic == 3) {
assign_val(prec_ewald, "PREC_EWALD", general_var_dict);
if (prec_ewald <= 0.0 || prec_ewald >= 1.0) {
Expand Down Expand Up @@ -314,6 +320,7 @@ void Input::parse_general_vars()
dynamical->na_sigma = na_sigma;
writes->nbands = nbands;
dynamical->file_born = borninfo;
dynamical->band_connection = band_connection;
integration->epsilon = epsilon;
fcs_phonon->file_fcs = fcsinfo;
if (!fc2info.empty()) {
Expand Down
47 changes: 42 additions & 5 deletions anphon/write_phonons.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,7 @@ void Writes::write_input_vars()
<< "; EPSILON = " << integration->epsilon << std::endl;
std::cout << std::endl;
std::cout << " CLASSICAL = " << thermodynamics->classical << std::endl;
std::cout << " BCONNECT = " << dynamical->band_connection << std::endl;
std::cout << std::endl;

if (phon->mode == "RTA") {
Expand Down Expand Up @@ -597,18 +598,54 @@ void Writes::write_phonon_bands()
ofs_bands << "#" << str_kval << std::endl;
ofs_bands << "# k-axis, Eigenvalues [cm^-1]" << std::endl;

for (i = 0; i < nk; ++i) {
ofs_bands << std::setw(8) << std::fixed << kaxis[i];
for (j = 0; j < nbands; ++j) {
ofs_bands << std::setw(15) << std::scientific << in_kayser(eval[i][j]);
if (dynamical->band_connection == 0) {
for (i = 0; i < nk; ++i) {
ofs_bands << std::setw(8) << std::fixed << kaxis[i];
for (j = 0; j < nbands; ++j) {
ofs_bands << std::setw(15) << std::scientific << in_kayser(eval[i][j]);
}
ofs_bands << std::endl;
}
} else {
for (i = 0; i < nk; ++i) {
ofs_bands << std::setw(8) << std::fixed << kaxis[i];
for (j = 0; j < nbands; ++j) {
ofs_bands << std::setw(15) << std::scientific
<< in_kayser(eval[i][dynamical->index_bconnect[i][j]]);
}
ofs_bands << std::endl;
}
ofs_bands << std::endl;
}

ofs_bands.close();

std::cout << " " << std::setw(input->job_title.length() + 12) << std::left << file_bands;
std::cout << " : Phonon band structure" << std::endl;

if (dynamical->band_connection == 2) {
std::ofstream ofs_connect;
std::string file_connect = input->job_title + ".connection";

ofs_connect.open(file_connect.c_str(), std::ios::out);
if (!ofs_connect)
error->exit("write_phonon_bands",
"cannot open file_connect");

ofs_connect << "# " << str_kpath << std::endl;
ofs_connect << "#" << str_kval << std::endl;
ofs_connect << "# k-axis, mapping" << std::endl;

for (i = 0; i < nk; ++i) {
ofs_connect << std::setw(8) << std::fixed << kaxis[i];
for (j = 0; j < nbands; ++j) {
ofs_connect << std::setw(5) << dynamical->index_bconnect[i][j] + 1;
}
ofs_connect << std::endl;
}
ofs_connect.close();
std::cout << " " << std::setw(input->job_title.length() + 12) << std::left << file_connect;
std::cout << " : Connectivity map information of band dispersion" << std::endl;
}
}

void Writes::write_phonon_vel()
Expand Down

0 comments on commit 3f41331

Please sign in to comment.