diff --git a/alm/alamode.cpp b/alm/alamode.cpp index 46244c24..5f1a6616 100644 --- a/alm/alamode.cpp +++ b/alm/alamode.cpp @@ -57,6 +57,7 @@ ALM::ALM(int narg, char **arg) if (mode == "fitting") { + fcs->init(); constraint->setup(); fitting->fitmain(); writes->writeall(); @@ -94,7 +95,6 @@ void ALM::initialize() files->init(); symmetry->init(); interaction->init(); - fcs->init(); } ALM::~ALM() diff --git a/alm/alamode.h b/alm/alamode.h index 3a2771c3..7fdbb494 100644 --- a/alm/alamode.h +++ b/alm/alamode.h @@ -40,4 +40,3 @@ namespace ALM_NS std::string mode; }; } - diff --git a/alm/combination.h b/alm/combination.h index c0c8f93c..492f4c13 100644 --- a/alm/combination.h +++ b/alm/combination.h @@ -25,7 +25,9 @@ namespace ALM_NS public: - CombinationWithRepetition() {}; + CombinationWithRepetition() + { + }; template CombinationWithRepetition(InputIter begin, @@ -75,4 +77,3 @@ namespace ALM_NS } }; } - diff --git a/alm/constraint.cpp b/alm/constraint.cpp index 85195e5b..4f0afe68 100644 --- a/alm/constraint.cpp +++ b/alm/constraint.cpp @@ -21,6 +21,7 @@ #include "error.h" #include #include "mathfunctions.h" +#include #ifdef _USE_EIGEN #include @@ -30,22 +31,33 @@ using namespace ALM_NS; Constraint::Constraint(ALM *alm) : Pointers(alm) { + const_symmetry = nullptr; + const_fix = nullptr; + const_relate = nullptr; + index_bimap = nullptr; + const_mat = nullptr; + const_rhs = nullptr; } Constraint::~Constraint() { - if (exist_constraint && alm->mode == "fitting") { - + if (const_symmetry) { memory->deallocate(const_symmetry); - - if (constraint_algebraic) { - memory->deallocate(const_fix); - memory->deallocate(const_relate); - memory->deallocate(index_bimap); - } else { - memory->deallocate(const_mat); - memory->deallocate(const_rhs); - } + } + if (const_fix) { + memory->deallocate(const_fix); + } + if (const_relate) { + memory->deallocate(const_relate); + } + if (index_bimap) { + memory->deallocate(index_bimap); + } + if (const_mat) { + memory->deallocate(const_mat); + } + if (const_rhs) { + memory->deallocate(const_rhs); } } @@ -107,7 +119,7 @@ void Constraint::setup() extra_constraint_from_symmetry = false; memory->allocate(const_symmetry, interaction->maxorder); - constraint_from_symmetry(const_symmetry); + generate_symmetry_constraint_in_cartesian(const_symmetry); for (int order = 0; order < interaction->maxorder; ++order) { if (const_symmetry[order].size() > 0) extra_constraint_from_symmetry = true; } @@ -126,7 +138,7 @@ void Constraint::setup() int N = 0; for (i = 0; i < maxorder; ++i) { - N += fcs->ndup[i].size(); + N += fcs->nequiv[i].size(); } memory->allocate(const_translation, maxorder); @@ -181,7 +193,7 @@ void Constraint::setup() for (order = 0; order < maxorder; ++order) { - nparam = fcs->ndup[order].size(); + nparam = fcs->nequiv[order].size(); memory->allocate(arr_tmp, nparam); for (it_const = const_translation[order].begin(); @@ -243,7 +255,7 @@ void Constraint::setup() memory->allocate(const_relate, maxorder); memory->allocate(index_bimap, maxorder); - get_mapping_constraint(maxorder, const_self, const_fix, + get_mapping_constraint(maxorder, fcs->nequiv, const_self, const_fix, const_relate, index_bimap, false); for (order = 0; order < maxorder; ++order) { @@ -260,11 +272,11 @@ void Constraint::setup() } if (fix_harmonic) { Pmax -= const_self[0].size(); - Pmax += fcs->ndup[0].size(); + Pmax += fcs->nequiv[0].size(); } if (fix_cubic) { Pmax -= const_self[1].size(); - Pmax += fcs->ndup[1].size(); + Pmax += fcs->nequiv[1].size(); } memory->allocate(const_mat, Pmax, N); memory->allocate(const_rhs, Pmax); @@ -302,13 +314,12 @@ void Constraint::calc_constraint_matrix(const int N, int &P) // Intra-order constraints for (order = 0; order < maxorder; ++order) { - int nparam = fcs->ndup[order].size(); + int nparam = fcs->nequiv[order].size(); if ((order == 0 && !fix_harmonic) || (order == 1 && !fix_cubic) || order > 1) { for (i = 0; i < N; ++i) arr_tmp[i] = 0.0; - for (std::vector::iterator p = const_self[order].begin(); - p != const_self[order].end(); ++p) { + for (auto p = const_self[order].begin(); p != const_self[order].end(); ++p) { ConstraintClass const_now = *p; for (i = 0; i < nparam; ++i) { arr_tmp[nshift + i] = const_now.w_const[i]; @@ -324,9 +335,9 @@ void Constraint::calc_constraint_matrix(const int N, int &P) int nshift2 = 0; for (order = 0; order < maxorder; ++order) { if (order > 0) { - int nparam2 = fcs->ndup[order - 1].size() + fcs->ndup[order].size(); + int nparam2 = fcs->nequiv[order - 1].size() + fcs->nequiv[order].size(); for (i = 0; i < N; ++i) arr_tmp[i] = 0.0; - for (std::vector::iterator p = const_rotation_cross[order].begin(); + for (auto p = const_rotation_cross[order].begin(); p != const_rotation_cross[order].end(); ++p) { ConstraintClass const_now = *p; for (i = 0; i < nparam2; ++i) { @@ -334,7 +345,7 @@ void Constraint::calc_constraint_matrix(const int N, int &P) } const_total.push_back(ConstraintClass(N, arr_tmp)); } - nshift2 += fcs->ndup[order - 1].size(); + nshift2 += fcs->nequiv[order - 1].size(); } } memory->deallocate(arr_tmp); @@ -349,7 +360,7 @@ void Constraint::calc_constraint_matrix(const int N, int &P) std::cout << " of the reference " << fc2_file << std::endl; std::cout << " Constraint for HARMONIC IFCs will be updated accordingly." << std::endl << std::endl; - P += fcs->ndup[0].size(); + P += fcs->nequiv[0].size(); } if (fix_cubic) { @@ -357,7 +368,7 @@ void Constraint::calc_constraint_matrix(const int N, int &P) std::cout << " of the reference " << fc3_file << std::endl; std::cout << " Constraint for ANHARM3 IFCs will be updated accordingly." << std::endl << std::endl; - P += fcs->ndup[1].size(); + P += fcs->nequiv[1].size(); } for (i = 0; i < P; ++i) { @@ -372,7 +383,7 @@ void Constraint::calc_constraint_matrix(const int N, int &P) if (fix_harmonic) { double *const_rhs_tmp; - int nfcs_tmp = fcs->ndup[0].size(); + int nfcs_tmp = fcs->nequiv[0].size(); memory->allocate(const_rhs_tmp, nfcs_tmp); system->load_reference_system_xml(fc2_file, 0, const_rhs_tmp); @@ -387,9 +398,9 @@ void Constraint::calc_constraint_matrix(const int N, int &P) } if (fix_cubic) { - int ishift = fcs->ndup[0].size(); + int ishift = fcs->nequiv[0].size(); double *const_rhs_tmp; - int nfcs_tmp = fcs->ndup[1].size(); + int nfcs_tmp = fcs->nequiv[1].size(); memory->allocate(const_rhs_tmp, nfcs_tmp); system->load_reference_system_xml(fc3_file, 1, const_rhs_tmp); @@ -403,8 +414,7 @@ void Constraint::calc_constraint_matrix(const int N, int &P) memory->deallocate(const_rhs_tmp); } - for (std::vector::iterator p = const_total.begin(); - p != const_total.end(); ++p) { + for (auto p = const_total.begin(); p != const_total.end(); ++p) { for (i = 0; i < N; ++i) { const_mat[irow][i] = (*p).w_const[i]; } @@ -415,6 +425,7 @@ void Constraint::calc_constraint_matrix(const int N, int &P) void Constraint::get_mapping_constraint(const int nmax, + std::vector *nequiv, std::vector *const_in, std::vector *const_fix_out, std::vector *const_relate_out, @@ -423,7 +434,7 @@ void Constraint::get_mapping_constraint(const int nmax, { int order; unsigned int i; - double const_tol = eps8; + //double const_tol = eps8; bool *fix_forceconstant; std::string *file_forceconstant; @@ -446,7 +457,7 @@ void Constraint::get_mapping_constraint(const int nmax, int nparam; for (order = 0; order < nmax; ++order) { - nparam = fcs->ndup[order].size(); + nparam = nequiv[order].size(); if (fix_forceconstant[order]) { @@ -466,11 +477,10 @@ void Constraint::get_mapping_constraint(const int nmax, std::vector alpha_tmp; std::vector p_index_tmp; - for (std::vector::reverse_iterator p = const_in[order].rbegin(); - p != const_in[order].rend(); ++p) { + for (auto p = const_in[order].rbegin(); p != const_in[order].rend(); ++p) { p_index_target = -1; for (i = 0; i < nparam; ++i) { - if (std::abs((*p).w_const[i]) > const_tol) { + if (std::abs((*p).w_const[i]) > tolerance_constraint) { p_index_target = i; break; } @@ -485,7 +495,7 @@ void Constraint::get_mapping_constraint(const int nmax, p_index_tmp.clear(); for (i = p_index_target + 1; i < nparam; ++i) { - if (std::abs((*p).w_const[i]) > const_tol) { + if (std::abs((*p).w_const[i]) > tolerance_constraint) { alpha_tmp.push_back((*p).w_const[i]); p_index_tmp.push_back(i); } @@ -509,7 +519,7 @@ void Constraint::get_mapping_constraint(const int nmax, for (order = 0; order < nmax; ++order) { - nparam = fcs->ndup[order].size(); + nparam = nequiv[order].size(); for (i = 0; i < nparam; ++i) { has_constraint[order].push_back(0); @@ -527,7 +537,7 @@ void Constraint::get_mapping_constraint(const int nmax, int icount; for (order = 0; order < nmax; ++order) { - nparam = fcs->ndup[order].size(); + nparam = nequiv[order].size(); icount = 0; for (i = 0; i < nparam; ++i) { @@ -541,193 +551,248 @@ void Constraint::get_mapping_constraint(const int nmax, } memory->deallocate(has_constraint); + memory->deallocate(fix_forceconstant); + memory->deallocate(file_forceconstant); } - -void Constraint::constraint_from_symmetry(std::vector *const_out) +void Constraint::generate_symmetry_constraint_in_cartesian(std::vector *const_out) { // Create constraint matrices arising from the crystal symmetry. - int i; unsigned int isym; - int ixyz, nxyz; int order; int maxorder = interaction->maxorder; - - int *index_tmp; - int **xyzcomponent; - int nparams; - - double *arr_constraint; bool has_constraint_from_symm = false; - std::set list_found; - std::vector> const_mat; for (isym = 0; isym < symmetry->nsym; ++isym) { - if (symmetry->sym_available[isym]) continue; + if (symmetry->SymmData[isym].compatible_with_cartesian) continue; has_constraint_from_symm = true; } - for (order = 0; order < maxorder; ++order) const_out[order].clear(); - if (has_constraint_from_symm) { std::cout << " Generating constraints from crystal symmetry ..." << std::endl; } - memory->allocate(index_tmp, maxorder + 1); + for (order = 0; order < maxorder; ++order) { + if (has_constraint_from_symm) { + std::cout << " " << std::setw(8) << interaction->str_order[order] << " ..."; + } + get_symmetry_constraint(order, interaction->pairs[order], + symmetry->SymmData, "Cartesian", + fcs->fc_table[order], fcs->nequiv[order], + const_out[order]); + if (has_constraint_from_symm) { + std::cout << " done." << std::endl; + } + } + if (has_constraint_from_symm) { + std::cout << " Finished !" << std::endl << std::endl; + } +} + +void Constraint::get_symmetry_constraint(const int order, const std::set pairs, + const std::vector symmop, + const std::string basis, + const std::vector fc_table, + const std::vector nequiv, + std::vector &const_out) +{ + // Create constraint matrices arising from the crystal symmetry. + + int i, j; + unsigned int isym; + int ixyz, nxyz; + int *index_tmp; + int **xyzcomponent; + int nparams; + int counter; + int nsym_in_use; + double *arr_constraint; + bool has_constraint_from_symm = false; + std::unordered_set list_found; + std::vector> const_mat; + int **map_sym; + double ***rotation; + + if (order < 0) return; + + int nsym = symmop.size(); + int nat = system->nat; + nparams = nequiv.size(); + + if (nparams == 0) return; + + memory->allocate(rotation, nsym, 3, 3); + memory->allocate(map_sym, nat, nsym); + memory->allocate(index_tmp, order + 2); + nxyz = static_cast(std::pow(static_cast(3), order + 2)); + memory->allocate(xyzcomponent, nxyz, order + 2); + fcs->get_xyzcomponent(order + 2, xyzcomponent); + nsym_in_use = 0; + counter = 0; const_mat.clear(); + int nfcs = fc_table.size(); - for (order = 0; order < maxorder; ++order) { + if (basis == "Cartesian") { - nparams = fcs->ndup[order].size(); + for (auto it = symmop.begin(); it != symmop.end(); ++it) { - if (has_constraint_from_symm) { - std::cout << " " << std::setw(8) << interaction->str_order[order] << " ..."; - if (nparams == 0) { - std::cout << " No parameters! Skipped." << std::endl; - continue; + if (!(*it).compatible_with_cartesian) { + + for (i = 0; i < 3; ++i) { + for (j = 0; j < 3; ++j) { + rotation[nsym_in_use][i][j] = (*it).rotation_cart[i][j]; + } + } + for (i = 0; i < nat; ++i) { + map_sym[i][nsym_in_use] = symmetry->map_sym[i][counter]; + } + ++nsym_in_use; } + ++counter; } - // Generate temporary list of parameters - list_found.clear(); - for (auto p = fcs->fc_set[order].begin(); p != fcs->fc_set[order].end(); ++p) { - for (i = 0; i < order + 2; ++i) index_tmp[i] = (*p).elems[i]; - list_found.insert(FcProperty(order + 2, (*p).coef, - index_tmp, (*p).mother)); + } else if (basis == "Lattice") { + + for (auto it = symmop.begin(); it != symmop.end(); ++it) { + if (!(*it).compatible_with_lattice) { + for (i = 0; i < 3; ++i) { + for (j = 0; j < 3; ++j) { + rotation[nsym_in_use][i][j] + = static_cast((*it).rotation[i][j]); + } + } + for (i = 0; i < nat; ++i) { + map_sym[i][nsym_in_use] = symmetry->map_sym[i][counter]; + } + ++nsym_in_use; + } + ++counter; } - nxyz = static_cast(std::pow(static_cast(3), order + 2)); - memory->allocate(xyzcomponent, nxyz, order + 2); - fcs->get_xyzcomponent(order + 2, xyzcomponent); - int nfcs = fcs->fc_set[order].size(); + } else { + memory->deallocate(rotation); + memory->deallocate(map_sym); + error->exit("get_symmetry_constraint", "Invalid basis input"); + } + + // Generate temporary list of parameters + list_found.clear(); + for (auto p = fc_table.begin(); p != fc_table.end(); ++p) { + for (i = 0; i < order + 2; ++i) index_tmp[i] = (*p).elems[i]; + list_found.insert(FcProperty(order + 2, (*p).sign, + index_tmp, (*p).mother)); + } + #ifdef _OPENMP #pragma omp parallel #endif - { - int j; - int i_prim; - int loc_nonzero; - int *ind; - int *atm_index, *atm_index_symm; - int *xyz_index; - double c_tmp; - - std::set::iterator iter_found; - std::vector const_now_omp; - std::vector> const_omp; - - memory->allocate(ind, order + 2); - memory->allocate(atm_index, order + 2); - memory->allocate(atm_index_symm, order + 2); - memory->allocate(xyz_index, order + 2); - - const_omp.clear(); - const_now_omp.resize(nparams); + { + int j; + int i_prim; + int loc_nonzero; + int *ind; + int *atm_index, *atm_index_symm; + int *xyz_index; + double c_tmp; + + std::unordered_set::iterator iter_found; + std::vector const_now_omp; + std::vector> const_omp; + + memory->allocate(ind, order + 2); + memory->allocate(atm_index, order + 2); + memory->allocate(atm_index_symm, order + 2); + memory->allocate(xyz_index, order + 2); + + const_omp.clear(); + const_now_omp.resize(nparams); #ifdef _OPENMP -#pragma omp for private(i, isym, ixyz) +#pragma omp for private(i, isym, ixyz), schedule(static) #endif - for (int ii = 0; ii < nfcs; ++ii) { - FcProperty list_tmp = fcs->fc_set[order][ii]; - - for (i = 0; i < order + 2; ++i) { - atm_index[i] = list_tmp.elems[i] / 3; - xyz_index[i] = list_tmp.elems[i] % 3; - } + for (int ii = 0; ii < nfcs; ++ii) { + FcProperty list_tmp = fc_table[ii]; - for (isym = 0; isym < symmetry->nsym; ++isym) { - - if (symmetry->sym_available[isym]) continue; + for (i = 0; i < order + 2; ++i) { + atm_index[i] = list_tmp.elems[i] / 3; + xyz_index[i] = list_tmp.elems[i] % 3; + } - for (i = 0; i < order + 2; ++i) - atm_index_symm[i] = symmetry->map_sym[atm_index[i]][isym]; - if (!fcs->is_inprim(order + 2, atm_index_symm)) continue; + for (isym = 0; isym < nsym_in_use; ++isym) { - for (i = 0; i < nparams; ++i) const_now_omp[i] = 0.0; + for (i = 0; i < order + 2; ++i) + atm_index_symm[i] = map_sym[atm_index[i]][isym]; + if (!fcs->is_inprim(order + 2, atm_index_symm)) continue; - const_now_omp[list_tmp.mother] = -list_tmp.coef; + for (i = 0; i < nparams; ++i) const_now_omp[i] = 0.0; - for (ixyz = 0; ixyz < nxyz; ++ixyz) { - for (i = 0; i < order + 2; ++i) - ind[i] = 3 * atm_index_symm[i] + xyzcomponent[ixyz][i]; + const_now_omp[list_tmp.mother] = -list_tmp.sign; - i_prim = fcs->min_inprim(order + 2, ind); - std::swap(ind[0], ind[i_prim]); - fcs->sort_tail(order + 2, ind); + for (ixyz = 0; ixyz < nxyz; ++ixyz) { + for (i = 0; i < order + 2; ++i) + ind[i] = 3 * atm_index_symm[i] + xyzcomponent[ixyz][i]; - iter_found = list_found.find(FcProperty(order + 2, 1.0, ind, 1)); - if (iter_found != list_found.end()) { - c_tmp = fcs->coef_sym(order + 2, isym, xyz_index, xyzcomponent[ixyz]); - const_now_omp[(*iter_found).mother] += (*iter_found).coef * c_tmp; - } - } + i_prim = fcs->min_inprim(order + 2, ind); + std::swap(ind[0], ind[i_prim]); + fcs->sort_tail(order + 2, ind); - if (!is_allzero(const_now_omp, loc_nonzero)) { - if (const_now_omp[loc_nonzero] < 0.0) { - for (j = 0; j < nparams; ++j) const_now_omp[j] *= -1.0; - } - const_omp.push_back(const_now_omp); + iter_found = list_found.find(FcProperty(order + 2, 1.0, ind, 1)); + if (iter_found != list_found.end()) { + c_tmp = fcs->coef_sym(order + 2, rotation[isym], xyz_index, xyzcomponent[ixyz]); + const_now_omp[(*iter_found).mother] += (*iter_found).sign * c_tmp; } - - } // close isym loop - - - // sort-->uniq the array - std::sort(const_omp.begin(), const_omp.end()); - const_omp.erase(std::unique(const_omp.begin(), const_omp.end()), - const_omp.end()); - - // Merge vectors -#pragma omp critical - { - for (std::vector>::iterator it = const_omp.begin(); - it != const_omp.end(); ++it) { - const_mat.push_back(*it); + } + if (!is_allzero(const_now_omp, eps8, loc_nonzero)) { + if (const_now_omp[loc_nonzero] < 0.0) { + for (j = 0; j < nparams; ++j) const_now_omp[j] *= -1.0; } + const_omp.push_back(const_now_omp); } - const_omp.clear(); - } // close ii loop + } // close isym loop - memory->deallocate(ind); - memory->deallocate(atm_index); - memory->deallocate(atm_index_symm); - memory->deallocate(xyz_index); + if (const_omp.size() > nparams) rref(const_omp, tolerance_constraint); - } // close openmp region + } // close ii loop - memory->allocate(arr_constraint, nparams); - for (std::vector>::reverse_iterator it = const_mat.rbegin(); - it != const_mat.rend(); ++it) { - for (i = 0; i < (*it).size(); ++i) { - arr_constraint[i] = (*it)[i]; + memory->deallocate(ind); + memory->deallocate(atm_index); + memory->deallocate(atm_index_symm); + memory->deallocate(xyz_index); + +#pragma omp critical + { + for (auto it = const_omp.begin(); it != const_omp.end(); ++it) { + const_mat.push_back(*it); } - const_out[order].push_back(ConstraintClass(nparams, - arr_constraint)); } - const_mat.clear(); + const_omp.clear(); + } // close openmp region - memory->deallocate(xyzcomponent); - memory->deallocate(arr_constraint); - - remove_redundant_rows(nparams, const_out[order], eps8); - - if (has_constraint_from_symm) { - std::cout << " done." << std::endl; + memory->allocate(arr_constraint, nparams); + for (auto it = const_mat.crbegin(); it != const_mat.crend(); ++it) { + for (i = 0; i < nparams; ++i) { + arr_constraint[i] = (*it)[i]; } - } // close loop order + const_out.push_back(ConstraintClass(nparams, + arr_constraint)); + } + const_mat.clear(); + memory->deallocate(xyzcomponent); + memory->deallocate(arr_constraint); memory->deallocate(index_tmp); + memory->deallocate(rotation); + memory->deallocate(map_sym); - if (has_constraint_from_symm) { - std::cout << " Finished !" << std::endl << std::endl; - } + remove_redundant_rows(nparams, const_out, tolerance_constraint); } + void Constraint::translational_invariance() { // Create constraint matrix for the translational invariance (aka acoustic sum rule). @@ -744,7 +809,7 @@ void Constraint::translational_invariance() int **xyzcomponent; int ixyz, nxyz; - int natmin = symmetry->natmin; + int natmin = symmetry->nat_prim; int nat = system->nat; int nparams; @@ -752,8 +817,8 @@ void Constraint::translational_invariance() double *arr_constraint; std::vector intlist, data; - std::set list_found; - std::set::iterator iter_found; + std::unordered_set list_found; + std::unordered_set::iterator iter_found; std::vector> data_vec; std::vector list_vec; std::vector::iterator iter_vec; @@ -770,7 +835,7 @@ void Constraint::translational_invariance() const_mat.clear(); - nparams = fcs->ndup[order].size(); + nparams = fcs->nequiv[order].size(); if (nparams == 0) { std::cout << " No parameters! Skipped." << std::endl; @@ -781,16 +846,15 @@ void Constraint::translational_invariance() list_found.clear(); - for (std::vector::iterator p = fcs->fc_set[order].begin(); - p != fcs->fc_set[order].end(); ++p) { + for (auto p = fcs->fc_table[order].cbegin(); p != fcs->fc_table[order].cend(); ++p) { for (i = 0; i < order + 2; ++i) { ind[i] = (*p).elems[i]; } - if (list_found.find(FcProperty(order + 2, (*p).coef, + if (list_found.find(FcProperty(order + 2, (*p).sign, ind, (*p).mother)) != list_found.end()) { error->exit("translational invariance", "Duplicate interaction list found"); } - list_found.insert(FcProperty(order + 2, (*p).coef, + list_found.insert(FcProperty(order + 2, (*p).sign, ind, (*p).mother)); } @@ -828,10 +892,10 @@ void Constraint::translational_invariance() iter_found = list_found.find( FcProperty(order + 2, 1.0, intarr, 1)); - // If found a IFC + // If found an IFC if (iter_found != list_found.end()) { // Round the coefficient to integer - const_now[(*iter_found).mother] += nint((*iter_found).coef); + const_now[(*iter_found).mother] += nint((*iter_found).sign); } } @@ -932,7 +996,7 @@ void Constraint::translational_invariance() iter_found = list_found.find(FcProperty(order + 2, 1.0, intarr_copy_omp, 1)); if (iter_found != list_found.end()) { - const_now_omp[(*iter_found).mother] += nint((*iter_found).coef); + const_now_omp[(*iter_found).mother] += nint((*iter_found).sign); } } @@ -955,8 +1019,7 @@ void Constraint::translational_invariance() // Merge vectors #pragma omp critical { - for (std::vector>::iterator it = const_omp.begin(); - it != const_omp.end(); ++it) { + for (auto it = const_omp.begin(); it != const_omp.end(); ++it) { const_mat.push_back(*it); } } @@ -982,11 +1045,10 @@ void Constraint::translational_invariance() memory->deallocate(xyzcomponent); memory->deallocate(intarr); memory->deallocate(intarr_copy); - // Copy to constraint class + // Copy to constraint class const_translation[order].clear(); - for (std::vector>::reverse_iterator it = const_mat.rbegin(); - it != const_mat.rend(); ++it) { + for (auto it = const_mat.rbegin(); it != const_mat.rend(); ++it) { for (i = 0; i < (*it).size(); ++i) { arr_constraint[i] = static_cast((*it)[i]); } @@ -1021,7 +1083,7 @@ void Constraint::rotational_invariance() int icrd, jcrd; int order; int maxorder = interaction->maxorder; - int natmin = symmetry->natmin; + int natmin = symmetry->nat_prim; int mu, nu; int ixyz, nxyz, nxyz2; int mu_lambda, lambda; @@ -1042,9 +1104,9 @@ void Constraint::rotational_invariance() std::vector interaction_list, interaction_list_old, interaction_list_now; - std::set list_found; - std::set list_found_last; - std::set::iterator iter_found; + std::unordered_set list_found; + std::unordered_set list_found_last; + std::unordered_set::iterator iter_found; CombinationWithRepetition g; @@ -1059,7 +1121,7 @@ void Constraint::rotational_invariance() for (order = 0; order < maxorder; ++order) { - nparams[order] = fcs->ndup[order].size(); + nparams[order] = fcs->nequiv[order].size(); if (order == 0) { std::cout << " Constraints between " << std::setw(8) @@ -1088,12 +1150,11 @@ void Constraint::rotational_invariance() list_found.clear(); - for (std::vector::iterator p = fcs->fc_set[order].begin(); - p != fcs->fc_set[order].end(); ++p) { + for (auto p = fcs->fc_table[order].begin(); p != fcs->fc_table[order].end(); ++p) { for (i = 0; i < order + 2; ++i) { ind[i] = (*p).elems[i]; } - list_found.insert(FcProperty(order + 2, (*p).coef, + list_found.insert(FcProperty(order + 2, (*p).sign, ind, (*p).mother)); } @@ -1127,7 +1188,7 @@ void Constraint::rotational_invariance() for (j = 0; j < nparam_sub; ++j) arr_constraint[j] = 0.0; - for (std::vector::iterator iter_list = interaction_list_now.begin(); + for (auto iter_list = interaction_list_now.begin(); iter_list != interaction_list_now.end(); ++iter_list) { jat = *iter_list; @@ -1163,8 +1224,7 @@ void Constraint::rotational_invariance() if (iter_found != list_found.end()) { - arr_constraint[(*iter_found).mother] - += (*iter_found).coef * vec_for_rot[nu]; + arr_constraint[(*iter_found).mother] += (*iter_found).sign * vec_for_rot[nu]; } // Exchange mu <--> nu and repeat again. @@ -1175,7 +1235,7 @@ void Constraint::rotational_invariance() interaction_index, 1)); if (iter_found != list_found.end()) { arr_constraint[(*iter_found).mother] - -= (*iter_found).coef * vec_for_rot[mu]; + -= (*iter_found).sign * vec_for_rot[mu]; } } @@ -1252,7 +1312,7 @@ void Constraint::rotational_invariance() for (j = 0; j < nparam_sub; ++j) arr_constraint[j] = 0.0; // Loop for m_{N+1}, a_{N+1} - for (std::vector::iterator iter_list = interaction_list.begin(); + for (auto iter_list = interaction_list.begin(); iter_list != interaction_list.end(); ++iter_list) { jat = *iter_list; @@ -1309,7 +1369,7 @@ void Constraint::rotational_invariance() iter_found = list_found.find(FcProperty(order + 2, 1.0, interaction_tmp, 1)); if (iter_found != list_found.end()) { arr_constraint[nparams[order - 1] + (*iter_found).mother] - += (*iter_found).coef * vec_for_rot[nu]; + += (*iter_found).sign * vec_for_rot[nu]; } // Exchange mu <--> nu and repeat again. @@ -1322,7 +1382,7 @@ void Constraint::rotational_invariance() iter_found = list_found.find(FcProperty(order + 2, 1.0, interaction_tmp, 1)); if (iter_found != list_found.end()) { arr_constraint[nparams[order - 1] + (*iter_found).mother] - -= (*iter_found).coef * vec_for_rot[mu]; + -= (*iter_found).sign * vec_for_rot[mu]; } } @@ -1350,7 +1410,7 @@ void Constraint::rotational_invariance() interaction_tmp, 1)); if (iter_found != list_found_last.end()) { arr_constraint[(*iter_found).mother] - += (*iter_found).coef * static_cast(levi_factor); + += (*iter_found).sign * static_cast(levi_factor); } } } @@ -1442,7 +1502,7 @@ void Constraint::rotational_invariance() interaction_tmp, 1)); if (iter_found != list_found.end()) { arr_constraint_self[(*iter_found).mother] - += (*iter_found).coef * static_cast(levi_factor); + += (*iter_found).sign * static_cast(levi_factor); } } // jcrd } // lambda @@ -1493,72 +1553,6 @@ void Constraint::remove_redundant_rows(const int n, std::vector &Constraint_vec, const double tolerance) { -#ifdef _USE_EIGEN_DISABLED - - // This function doesn't make the reduced row echelon form of the constraint matrix. - // It just returns the image of the matrix, though they are similar. - - using namespace Eigen; - int i; - int nrow = n; - int ncol = Constraint_vec.size(); - double *arr_tmp; - double **mat; - - if (ncol > 0) { - memory->allocate(arr_tmp, nrow); - MatrixXd mat_tmp(nrow, ncol); - - int icol = 0; - - for (std::vector::iterator p = Constraint_vec.begin(); p != Constraint_vec.end(); ++p) { - ConstraintClass const_now = *p; - for (i = 0; i < nrow; ++i) { - mat_tmp(i, icol) = const_now.w_const[i]; - } - ++icol; - } - - FullPivLU lu_decomp(mat_tmp); - lu_decomp.setThreshold(tolerance); - int nrank = lu_decomp.rank(); - MatrixXd c_reduced = lu_decomp.image(mat_tmp); - - memory->allocate(mat, nrank, nrow); - - i = 0; - for (icol = 0; icol < nrank; ++icol) { - for (int irow = 0; irow < nrow; ++irow) { - mat[icol][irow] = c_reduced(irow, icol); - } - } - // std::cout << "nrank = " << nrank << std::flush << std::endl; - rref(nrank, nrow, mat, nrank, tolerance); - - Constraint_vec.clear(); - - for (i = 0; i < nrank; ++i) { - for (int j = 0; j < i; ++j) arr_tmp[j] = 0.0; - - for (int j = i; j < nrow; ++j) { - arr_tmp[j] = mat[i][j]; - } - Constraint_vec.push_back(ConstraintClass(nrow, arr_tmp)); - } - - memory->deallocate(mat); - - // for (icol = 0; icol < nrank; ++icol) { - // for (int irow = 0; irow < nrow; ++irow) { - // arr_tmp[irow] = c_reduced(irow, icol); - // } - - // Constraint_Set.insert(ConstraintClass(nrow, arr_tmp)); - // } - - memory->deallocate(arr_tmp); - } -#else int i, j; int nparam = n; @@ -1574,8 +1568,7 @@ void Constraint::remove_redundant_rows(const int n, i = 0; - for (std::vector::iterator p = Constraint_vec.begin(); - p != Constraint_vec.end(); ++p) { + for (auto p = Constraint_vec.begin(); p != Constraint_vec.end(); ++p) { for (j = 0; j < nparam; ++j) { mat_tmp[i][j] = (*p).w_const[j]; } @@ -1584,38 +1577,6 @@ void Constraint::remove_redundant_rows(const int n, rref(nconst, nparam, mat_tmp, nrank, tolerance); - /* - // // Transpose matrix A - - memory->allocate(arr_tmp, nconst * nparam); - - k = 0; - - for (j = 0; j < nparam; ++j) { - for (i = 0; i < nconst; ++i) { - arr_tmp[k++] = mat_tmp[i][j]; - } - } - - // Perform LU decomposition - - int nmin = std::min(nconst, nparam); - memory->allocate(ipiv, nmin); - - dgetrf_(&nconst, &nparam, arr_tmp, &nconst, ipiv, &INFO); - - k = 0; - - for (j = 0; j < nparam; ++j) { - for (i = 0; i < nconst; ++i) { - mat_tmp[i][j] = arr_tmp[k++]; - } - } - - memory->deallocate(arr_tmp); - memory->deallocate(ipiv); - */ - memory->allocate(arr_tmp, nparam); Constraint_vec.clear(); @@ -1624,18 +1585,22 @@ void Constraint::remove_redundant_rows(const int n, for (j = 0; j < i; ++j) arr_tmp[j] = 0.0; for (j = i; j < nparam; ++j) { - arr_tmp[j] = mat_tmp[i][j]; + if (std::abs(mat_tmp[i][j]) < tolerance) { + arr_tmp[j] = 0.0; + } else { + arr_tmp[j] = mat_tmp[i][j]; + } } Constraint_vec.push_back(ConstraintClass(nparam, arr_tmp)); } memory->deallocate(mat_tmp); memory->deallocate(arr_tmp); - } -#endif + } } + int Constraint::levi_civita(const int i, const int j, const int k) { int epsilon = (j - i) * (k - i) * (k - j) / 2; @@ -1664,11 +1629,11 @@ bool Constraint::is_allzero(const std::vector vec, int &loc) return true; } -bool Constraint::is_allzero(const std::vector vec, int &loc) +bool Constraint::is_allzero(const std::vector vec, const double tol, int &loc) { loc = -1; for (int i = 0; i < vec.size(); ++i) { - if (std::abs(vec[i]) > eps) { + if (std::abs(vec[i]) > tol) { loc = i; return false; } @@ -1735,9 +1700,7 @@ void Constraint::rref(int nrows, int irow, icol, jrow, jcol; int pivot; - double tmp, *arr; - - memory->allocate(arr, ncols); + double tmp; nrank = 0; @@ -1763,7 +1726,7 @@ void Constraint::rref(int nrows, if (std::abs(mat[pivot][icol]) > tolerance) ++nrank; if (pivot != irow) { -#pragma omp parallel for private(tmp) + //#pragma omp parallel for private(tmp) for (jcol = icol; jcol < ncols; ++jcol) { tmp = mat[pivot][jcol]; mat[pivot][jcol] = mat[irow][jcol]; @@ -1773,7 +1736,7 @@ void Constraint::rref(int nrows, tmp = mat[irow][icol]; tmp = 1.0 / tmp; -#pragma omp parallel for + //#pragma omp parallel for for (jcol = icol; jcol < ncols; ++jcol) { mat[irow][jcol] *= tmp; } @@ -1782,12 +1745,183 @@ void Constraint::rref(int nrows, if (jrow == irow) continue; tmp = mat[jrow][icol]; -#pragma omp parallel for + //#pragma omp parallel for for (jcol = icol; jcol < ncols; ++jcol) { mat[jrow][jcol] -= tmp * mat[irow][jcol]; } } } +} + + +void Constraint::rref(std::vector> &mat, const double tolerance) +{ + // Return the reduced row echelon form (rref) of matrix mat. + // In addition, rank of the matrix is estimated. + + int irow, icol, jrow, jcol; + int pivot; + double tmp; + + int nrank = 0; + + icol = 0; + + int nrows = mat.size(); + + if (nrows == 0) return; + + int ncols = mat[0].size(); + + for (irow = 0; irow < nrows; ++irow) { + + pivot = irow; + + while (std::abs(mat[pivot][icol]) < tolerance) { + ++pivot; + + if (pivot == nrows) { + pivot = irow; + ++icol; + + if (icol == ncols) break; + } + } + + if (icol == ncols) break; + + if (std::abs(mat[pivot][icol]) > tolerance) ++nrank; + + if (pivot != irow) { + for (jcol = icol; jcol < ncols; ++jcol) { + tmp = mat[pivot][jcol]; + mat[pivot][jcol] = mat[irow][jcol]; + mat[irow][jcol] = tmp; + } + } + + tmp = mat[irow][icol]; + tmp = 1.0 / tmp; + for (jcol = icol; jcol < ncols; ++jcol) { + mat[irow][jcol] *= tmp; + } + + for (jrow = 0; jrow < nrows; ++jrow) { + if (jrow == irow) continue; + + tmp = mat[jrow][icol]; + for (jcol = icol; jcol < ncols; ++jcol) { + mat[jrow][jcol] -= tmp * mat[irow][jcol]; + } + } + } + + mat.erase(mat.begin() + nrank, mat.end()); + mat.shrink_to_fit(); +} + + +/* +void Constraint::rref_nofraction(std::vector> &mat) +{ + // Return the reduced row echelon form (rref) of matrix mat. + // In addition, rank of the matrix is estimated. + + int irow, icol, jrow, jcol; + int pivot; + int tmp, tmp2; + + int nrank = 0; + + icol = 0; + + int nrows = mat.size(); + int ncols = mat[0].size(); + + for (irow = 0; irow < nrows; ++irow) { + + pivot = irow; + + while (mat[pivot][icol] == 0) { + ++pivot; + + if (pivot == nrows) { + pivot = irow; + ++icol; + + if (icol == ncols) break; + } + } + + if (icol == ncols) break; + + if (std::abs(mat[pivot][icol]) > 0) ++nrank; + + // swap rows + if (pivot != irow) { +#pragma omp parallel for private(tmp) + for (jcol = icol; jcol < ncols; ++jcol) { + tmp = mat[pivot][jcol]; + mat[pivot][jcol] = mat[irow][jcol]; + mat[irow][jcol] = tmp; + } + } + + tmp = mat[irow][icol]; + + for (jrow = 0; jrow < nrows; ++jrow) { + if (jrow == irow) continue; + + tmp2 = mat[jrow][icol]; +#pragma omp parallel for + for (jcol = icol; jcol < ncols; ++jcol) { + mat[jrow][jcol] = mat[jrow][jcol] * tmp - tmp2 * mat[irow][jcol]; + } + } + } + + mat.erase(mat.begin() + nrank, mat.end()); + mat.shrink_to_fit(); +} +*/ + + +#ifdef _USE_EIGEN +void Constraint::get_column_space(std::vector> &mat) +{ +// Return the column space of matrix mat. + + using namespace Eigen; + + int irow, icol, jrow, jcol; + int pivot; + int tmp, tmp2; + + int nrank = 0; + + icol = 0; + + int nrows = mat.size(); + int ncols = mat[0].size(); + + MatrixXf A(ncols, nrows); + + for (irow = 0; irow < nrows; ++irow) { + for (icol = 0; icol < ncols; ++icol) { + A(icol, irow) = static_cast(mat[irow][icol]); + } + } + FullPivLU lu_decomp(A); + nrank = lu_decomp.rank(); + MatrixXf B = lu_decomp.image(A).transpose(); + + for (irow = 0; irow < nrank; ++irow) { + for (icol = 0; icol < ncols; ++icol) { + mat[irow][icol] = static_cast(B(irow, icol) + 0.5); + } + } - memory->deallocate(arr); + mat.erase(mat.begin() + nrank, mat.end()); + mat.shrink_to_fit(); } +#endif diff --git a/alm/constraint.h b/alm/constraint.h index c5c36015..ba5c3b27 100644 --- a/alm/constraint.h +++ b/alm/constraint.h @@ -16,6 +16,9 @@ #include #include "pointers.h" #include "constants.h" +#include "interaction.h" +#include "symmetry.h" +#include "fcs.h" #include namespace ALM_NS @@ -29,8 +32,8 @@ namespace ALM_NS ConstraintClass(const ConstraintClass &a) { - for (std::vector::const_iterator p = a.w_const.begin(); - p != a.w_const.end(); ++p) { + for (std::vector::const_iterator p = a.w_const.begin(); + p != a.w_const.end(); ++p) { w_const.push_back(*p); } } @@ -79,6 +82,19 @@ namespace ALM_NS } }; + inline bool equal_within_eps12(const std::vector &a, const std::vector &b) + { + int n = a.size(); + int m = b.size(); + if (n != m) return false; + double res = 0.0; + for (int i = 0; i < n; ++i) { + if (std::abs(a[i] - b[i]) > eps12) return false; + } + // if (std::sqrt(res)>eps12) return false; + return true; + } + class Constraint: protected Pointers { public: @@ -95,6 +111,7 @@ namespace ALM_NS double **const_mat; double *const_rhs; + double tolerance_constraint; bool exist_constraint; bool extra_constraint_from_symmetry; @@ -105,9 +122,16 @@ namespace ALM_NS std::vector *const_relate; boost::bimap *index_bimap; - void constraint_from_symmetry(std::vector *); + //void constraint_from_symmetry(std::vector *); + void get_symmetry_constraint(const int, const std::set, + const std::vector, + const std::string, + const std::vector, + const std::vector, + std::vector &); - void get_mapping_constraint(const int, std::vector *, + void get_mapping_constraint(const int, std::vector *, + std::vector *, std::vector *, std::vector *, boost::bimap *, const bool); @@ -131,19 +155,25 @@ namespace ALM_NS void setup_rotation_axis(bool [3][3]); bool is_allzero(const int, const double *, const int nshift = 0); bool is_allzero(const std::vector, int &); - bool is_allzero(const std::vector, int &); + bool is_allzero(const std::vector, const double, int &); void remove_redundant_rows(const int, std::vector &, const double tolerance = eps12); + // void remove_redundant_rows_integer(const int, std::vector> &); - void remove_redundant_rows2(const int, std::vector &, - const double tolerance = eps12); void rref(int, int, double **, int &, double tolerance = eps12); + void rref(std::vector> &, const double tolerance = eps12); + void rref_nofraction(std::vector> &); +#ifdef _USE_EIGEN + void get_column_space(std::vector> &); +#endif + + void generate_symmetry_constraint_in_cartesian(std::vector *); }; extern "C" { void dgetrf_(int *m, int *n, double *a, int *lda, int *ipiv, int *info); + void sgetrf_(int *m, int *n, float *a, int *lda, int *ipiv, int *info); } } - diff --git a/alm/error.h b/alm/error.h index 1b76c3bc..7d90525a 100644 --- a/alm/error.h +++ b/alm/error.h @@ -27,4 +27,3 @@ namespace ALM_NS void exit(const char *, const char *, const char *); }; } - diff --git a/alm/fcs.cpp b/alm/fcs.cpp index b49fa6cb..a2b4f54c 100644 --- a/alm/fcs.cpp +++ b/alm/fcs.cpp @@ -32,6 +32,11 @@ Fcs::Fcs(ALM *alm) : Pointers(alm) Fcs::~Fcs() { + if (alm->mode == "fitting") { + memory->deallocate(fc_table); + memory->deallocate(nequiv); + memory->deallocate(fc_zeros); + } }; void Fcs::init() @@ -42,53 +47,60 @@ void Fcs::init() std::cout << " FORCE CONSTANT" << std::endl; std::cout << " ==============" << std::endl << std::endl; - memory->allocate(nints, maxorder); - memory->allocate(nzero, maxorder); - memory->allocate(fc_set, maxorder); - memory->allocate(ndup, maxorder); + memory->allocate(fc_table, maxorder); + memory->allocate(nequiv, maxorder); + memory->allocate(fc_zeros, maxorder); - for (i = 0; i < maxorder; ++i) nzero[i] = 0; - generate_fclists(maxorder); + for (i = 0; i < maxorder; ++i) { + generate_force_constant_table(i, interaction->pairs[i], + symmetry->SymmData, "Cartesian", + fc_table[i], nequiv[i], fc_zeros[i], true); + } std::cout << std::endl; for (i = 0; i < maxorder; ++i) { std::cout << " Number of " << std::setw(9) << interaction->str_order[i] - << " FCs : " << ndup[i].size(); + << " FCs : " << nequiv[i].size(); std::cout << std::endl; } std::cout << std::endl; - // sort fc_set - - for (int order = 0; order < maxorder; ++order) { - if (ndup[order].size() > 0) { - std::sort(fc_set[order].begin(), - fc_set[order].begin() + ndup[order][0]); - int nbegin = ndup[order][0]; - int nend; - for (int mm = 1; mm < ndup[order].size(); ++mm) { - nend = nbegin + ndup[order][mm]; - std::sort(fc_set[order].begin() + nbegin, - fc_set[order].begin() + nend); - nbegin += ndup[order][mm]; - } - } - } - memory->deallocate(nints); - memory->deallocate(nzero); + // + // std::vector fc_test, fc_zeros_tmp; + // std::vector nmulti; + // generate_force_constant_table(1, interaction->pairs[1], + // symmetry->SymmData, "Lattice", + // fc_test, nmulti, fc_zeros_tmp, true); + // std::cout << "Nonzero independent IFCs:" << std::setw(5) << nmulti.size() << std::endl; + // std::cout << "Zero IFCs:" << std::endl; + // for (auto it = fc_zeros_tmp.begin(); it != fc_zeros_tmp.end(); ++it) { + // for (i = 0; i < (*it).elems.size(); ++i) { + // std::cout << std::setw(4) << (*it).elems[i]; + // } + // std::cout << std::setw(5) << (*it).mother; + // std::cout << std::endl; + // } + // std::cout << std::endl; + + timer->print_elapsed(); std::cout << " -------------------------------------------------------------------" << std::endl; std::cout << std::endl; } - -void Fcs::generate_fclists(int maxorder) +void Fcs::generate_force_constant_table(const int order, + const std::set pairs, + const std::vector symmop, + std::string basis, + std::vector &fc_vec, + std::vector &ndup, + std::vector &fc_zeros, + const bool store_zeros) { int i, j; int i1, i2; - int order; int i_prim; int *atmn, *atmn_mapped; int *ind, *ind_mapped; @@ -102,138 +114,184 @@ void Fcs::generate_fclists(int maxorder) int nmother; int nat = system->nat; + int nsym = symmop.size(); + int nsym_in_use; bool is_zero; bool *is_searched; + int counter; + int **map_sym; + double ***rotation; + + if (order < 0) return; + + memory->allocate(rotation, nsym, 3, 3); + memory->allocate(map_sym, nat, nsym); + nsym_in_use = 0; + counter = 0; + if (basis == "Cartesian") { + + for (auto it = symmop.begin(); it != symmop.end(); ++it) { + if ((*it).compatible_with_cartesian) { + for (i = 0; i < 3; ++i) { + for (j = 0; j < 3; ++j) { + rotation[nsym_in_use][i][j] = (*it).rotation_cart[i][j]; + } + } + for (i = 0; i < nat; ++i) { + map_sym[i][nsym_in_use] = symmetry->map_sym[i][counter]; + } + ++nsym_in_use; + } + ++counter; + } - std::cout << " Finding symmetrically-independent force constants ..." << std::endl; + } else if (basis == "Lattice") { - memory->allocate(atmn, maxorder + 1); - memory->allocate(atmn_mapped, maxorder + 1); - memory->allocate(ind, maxorder + 1); - memory->allocate(ind_mapped, maxorder + 1); - memory->allocate(ind_tmp, maxorder - 1); - memory->allocate(ind_mapped_tmp, maxorder + 1); - memory->allocate(is_searched, 3 * nat); + for (auto it = symmop.begin(); it != symmop.end(); ++it) { + if ((*it).compatible_with_lattice) { + for (i = 0; i < 3; ++i) { + for (j = 0; j < 3; ++j) { + rotation[nsym_in_use][i][j] + = static_cast((*it).rotation[i][j]); + } + } + for (i = 0; i < nat; ++i) { + map_sym[i][nsym_in_use] = symmetry->map_sym[i][counter]; + } + ++nsym_in_use; + } + ++counter; + } - for (order = 0; order < maxorder; ++order) { - std::cout << " " << std::setw(8) << interaction->str_order[order] << " ..."; + } else { + memory->deallocate(rotation); + memory->deallocate(map_sym); + error->exit("generate_force_constant_table", "Invalid basis inpout"); + } - fc_set[order].clear(); - ndup[order].clear(); - nmother = 0; + memory->allocate(atmn, order + 2); + memory->allocate(atmn_mapped, order + 2); + memory->allocate(ind, order + 2); + memory->allocate(ind_mapped, order + 2); + memory->allocate(ind_tmp, order); + memory->allocate(ind_mapped_tmp, order + 2); + memory->allocate(is_searched, 3 * nat); - nxyz = static_cast(std::pow(3.0, order + 2)); + fc_vec.clear(); + ndup.clear(); + fc_zeros.clear(); + nmother = 0; - memory->allocate(xyzcomponent, nxyz, order + 2); - get_xyzcomponent(order + 2, xyzcomponent); + nxyz = static_cast(std::pow(3.0, order + 2)); - std::set list_found; + memory->allocate(xyzcomponent, nxyz, order + 2); + get_xyzcomponent(order + 2, xyzcomponent); - for (std::set::iterator iter = interaction->pairs[order].begin(); - iter != interaction->pairs[order].end(); ++iter) { + std::set list_found; - for (i = 0; i < order + 2; ++i) atmn[i] = (*iter).iarray[i]; + for (auto iter = pairs.begin(); iter != pairs.end(); ++iter) { - for (i1 = 0; i1 < nxyz; ++i1) { - for (i = 0; i < order + 2; ++i) - ind[i] = 3 * atmn[i] + xyzcomponent[i1][i]; + for (i = 0; i < order + 2; ++i) atmn[i] = (*iter).iarray[i]; - if (!is_ascending(order + 2, ind)) continue; + for (i1 = 0; i1 < nxyz; ++i1) { + for (i = 0; i < order + 2; ++i) ind[i] = 3 * atmn[i] + xyzcomponent[i1][i]; - i_prim = min_inprim(order + 2, ind); - std::swap(ind[0], ind[i_prim]); - sort_tail(order + 2, ind); + if (!is_ascending(order + 2, ind)) continue; - is_zero = false; + i_prim = min_inprim(order + 2, ind); + std::swap(ind[0], ind[i_prim]); + sort_tail(order + 2, ind); - if (list_found.find(IntList(order + 2, ind)) != list_found.end()) continue; // Already exits! + is_zero = false; - // Search symmetrically-dependent parameter set + if (list_found.find(IntList(order + 2, ind)) != list_found.end()) continue; // Already exits! - int ndeps = 0; + // Search symmetrically-dependent parameter set - for (isym = 0; isym < symmetry->nsym; ++isym) { + int ndeps = 0; - if (!symmetry->sym_available[isym]) continue; + for (isym = 0; isym < nsym_in_use; ++isym) { - for (i = 0; i < order + 2; ++i) - atmn_mapped[i] = symmetry->map_sym[atmn[i]][isym]; + for (i = 0; i < order + 2; ++i) + atmn_mapped[i] = map_sym[atmn[i]][isym]; - if (!is_inprim(order + 2, atmn_mapped)) continue; + if (!is_inprim(order + 2, atmn_mapped)) continue; - for (i2 = 0; i2 < nxyz; ++i2) { - c_tmp = coef_sym(order + 2, isym, xyzcomponent[i1], xyzcomponent[i2]); - if (std::abs(c_tmp) > eps12) { - for (i = 0; i < order + 2; ++i) - ind_mapped[i] = 3 * atmn_mapped[i] + xyzcomponent[i2][i]; + for (i2 = 0; i2 < nxyz; ++i2) { + c_tmp = coef_sym(order + 2, rotation[isym], xyzcomponent[i1], xyzcomponent[i2]); + if (std::abs(c_tmp) > eps12) { + for (i = 0; i < order + 2; ++i) + ind_mapped[i] = 3 * atmn_mapped[i] + xyzcomponent[i2][i]; - i_prim = min_inprim(order + 2, ind_mapped); - std::swap(ind_mapped[0], ind_mapped[i_prim]); - sort_tail(order + 2, ind_mapped); + i_prim = min_inprim(order + 2, ind_mapped); + std::swap(ind_mapped[0], ind_mapped[i_prim]); + sort_tail(order + 2, ind_mapped); - if (!is_zero) { - bool zeroflag = true; - for (i = 0; i < order + 2; ++i) { - zeroflag = zeroflag & (ind[i] == ind_mapped[i]); - } - zeroflag = zeroflag & (std::abs(c_tmp + 1.0) < eps8); - is_zero = zeroflag; + if (!is_zero) { + bool zeroflag = true; + for (i = 0; i < order + 2; ++i) { + zeroflag = zeroflag & (ind[i] == ind_mapped[i]); } + zeroflag = zeroflag & (std::abs(c_tmp + 1.0) < eps8); + is_zero = zeroflag; + } - // Add to found list (set) and fcset (vector) if the created is new one. + // Add to found list (set) and fcset (vector) if the created is new one. - if (list_found.find(IntList(order + 2, ind_mapped)) == list_found.end()) { - list_found.insert(IntList(order + 2, ind_mapped)); + if (list_found.find(IntList(order + 2, ind_mapped)) == list_found.end()) { + list_found.insert(IntList(order + 2, ind_mapped)); - fc_set[order].push_back(FcProperty(order + 2, c_tmp, - ind_mapped, nmother)); - ++ndeps; + fc_vec.push_back(FcProperty(order + 2, c_tmp, + ind_mapped, nmother)); + ++ndeps; - // Add equivalent interaction list (permutation) if there are two or more indices - // which belong to the primitive cell. - // This procedure is necessary for fitting. + // Add equivalent interaction list (permutation) if there are two or more indices + // which belong to the primitive cell. + // This procedure is necessary for fitting. - for (i = 0; i < 3 * nat; ++i) is_searched[i] = false; - is_searched[ind_mapped[0]] = true; - for (i = 1; i < order + 2; ++i) { - if ((!is_searched[ind_mapped[i]]) && is_inprim(ind_mapped[i])) { + for (i = 0; i < 3 * nat; ++i) is_searched[i] = false; + is_searched[ind_mapped[0]] = true; + for (i = 1; i < order + 2; ++i) { + if ((!is_searched[ind_mapped[i]]) && is_inprim(ind_mapped[i])) { - for (j = 0; j < order + 2; ++j) ind_mapped_tmp[j] = ind_mapped[j]; - std::swap(ind_mapped_tmp[0], ind_mapped_tmp[i]); - sort_tail(order + 2, ind_mapped_tmp); - fc_set[order].push_back(FcProperty(order + 2, c_tmp, - ind_mapped_tmp, nmother)); + for (j = 0; j < order + 2; ++j) ind_mapped_tmp[j] = ind_mapped[j]; + std::swap(ind_mapped_tmp[0], ind_mapped_tmp[i]); + sort_tail(order + 2, ind_mapped_tmp); + fc_vec.push_back(FcProperty(order + 2, c_tmp, + ind_mapped_tmp, nmother)); - ++ndeps; + ++ndeps; - is_searched[ind_mapped[i]] = true; - } + is_searched[ind_mapped[i]] = true; } + } - } } } - } // close symmetry loop - - if (is_zero) { - for (i = 0; i < ndeps; ++i) fc_set[order].pop_back(); - ++nzero[order]; - } else { - ndup[order].push_back(ndeps); - ++nmother; } + } // close symmetry loop - } // close xyz component loop - } // close atom number loop (iterator) + if (is_zero) { + if (store_zeros) { + for (auto it = fc_vec.rbegin(); it != fc_vec.rbegin() + ndeps; ++it) { + (*it).mother = -1; + fc_zeros.push_back(*it); + } + } + for (i = 0; i < ndeps; ++i) fc_vec.pop_back(); + } else { + ndup.push_back(ndeps); + ++nmother; + } - memory->deallocate(xyzcomponent); - list_found.clear(); - std::cout << " done. " << std::endl; - } //close order loop + } // close xyz component loop + } // close atom number loop (iterator) + memory->deallocate(xyzcomponent); + list_found.clear(); memory->deallocate(atmn); memory->deallocate(atmn_mapped); memory->deallocate(ind); @@ -241,10 +299,24 @@ void Fcs::generate_fclists(int maxorder) memory->deallocate(ind_tmp); memory->deallocate(ind_mapped_tmp); memory->deallocate(is_searched); - - std::cout << " Finished!" << std::endl; + memory->deallocate(rotation); + memory->deallocate(map_sym); + + // sort fc_vec + + if (ndup.size() > 0) { + std::sort(fc_vec.begin(), fc_vec.begin() + ndup[0]); + int nbegin = ndup[0]; + int nend; + for (int mm = 1; mm < ndup.size(); ++mm) { + nend = nbegin + ndup[mm]; + std::sort(fc_vec.begin() + nbegin, fc_vec.begin() + nend); + nbegin += ndup[mm]; + } + } } + double Fcs::coef_sym(const int n, const int symnum, const int *arr1, @@ -254,7 +326,21 @@ double Fcs::coef_sym(const int n, int i; for (i = 0; i < n; ++i) { - tmp *= symmetry->symrel[symnum][arr2[i]][arr1[i]]; + tmp *= symmetry->SymmData[symnum].rotation_cart[arr2[i]][arr1[i]]; + } + return tmp; +} + +double Fcs::coef_sym(const int n, + double **rot, + const int *arr1, + const int *arr2) +{ + double tmp = 1.0; + int i; + + for (i = 0; i < n; ++i) { + tmp *= rot[arr2[i]][arr1[i]]; } return tmp; } @@ -271,7 +357,7 @@ bool Fcs::is_ascending(const int n, const int *arr) int Fcs::min_inprim(const int n, const int *arr) { int i, j, atmnum; - int natmin = symmetry->natmin; + int natmin = symmetry->nat_prim; int minloc; int *ind; @@ -307,7 +393,7 @@ int Fcs::min_inprim(const int n, const int *arr) bool Fcs::is_inprim(const int n, const int *arr) { int i, j; - int natmin = symmetry->natmin; + int natmin = symmetry->nat_prim; for (i = 0; i < n; ++i) { for (j = 0; j < natmin; ++j) { @@ -320,7 +406,7 @@ bool Fcs::is_inprim(const int n, const int *arr) bool Fcs::is_inprim(const int n) { int i, atmn; - int natmin = symmetry->natmin; + int natmin = symmetry->nat_prim; atmn = n / 3; diff --git a/alm/fcs.h b/alm/fcs.h index edf22db3..2eb69870 100644 --- a/alm/fcs.h +++ b/alm/fcs.h @@ -1,7 +1,7 @@ /* fcs.h - Copyright (c) 2014, 2015, 2016 Terumasa Tadano + Copyright (c) 2014--2017 Terumasa Tadano This file is distributed under the terms of the MIT license. Please see the file 'LICENCE.txt' in the root directory @@ -14,31 +14,32 @@ #include #include #include +#include "symmetry.h" +#include "interaction.h" namespace ALM_NS { class FcProperty { public: - std::vector elems; - double coef; + std::vector elems; // flattened index of (iatom, icoordinate) in the supercell + double sign; // factor (+1 or -1) to convert the mother FC to the child int mother; FcProperty(); FcProperty(const FcProperty &obj) { - coef = obj.coef; + sign = obj.sign; mother = obj.mother; - for (std::vector::const_iterator it = obj.elems.begin(); - it != obj.elems.end(); ++it) { + for (auto it = obj.elems.begin(); it != obj.elems.end(); ++it) { elems.push_back(*it); } }; FcProperty(const int n, const double c, const int *arr, const int m) { - coef = c; + sign = c; mother = m; for (int i = 0; i < n; ++i) { elems.push_back(arr[i]); @@ -50,6 +51,26 @@ namespace ALM_NS return std::lexicographical_compare(elems.begin(), elems.end(), a.elems.begin(), a.elems.end()); } + + bool operator==(const FcProperty &a) const + { + int n = elems.size(); + int n_ = a.elems.size(); + if (n != n_) return false; + for (int i = 0; i < n; ++i) { + if (elems[i] != a.elems[i]) return false; + } + return true; + } + }; + + class ForceConstantTable + { + public: + double fc_value; + int multiplicity; + std::vector fclist; + ForceConstantTable(); }; class Fcs: protected Pointers @@ -60,10 +81,9 @@ namespace ALM_NS void init(); - int *nzero; - - std::vector *ndup; - std::vector *fc_set; + std::vector *nequiv; + std::vector *fc_table; + std::vector *fc_zeros; std::string easyvizint(const int); void get_xyzcomponent(int, int **); @@ -73,11 +93,38 @@ namespace ALM_NS bool is_inprim(const int); int min_inprim(const int, const int *); double coef_sym(const int, const int, const int *, const int *); + double coef_sym(const int, double **, const int *, const int *); + + void generate_force_constant_table(const int, + const std::set, + const std::vector, + std::string, + std::vector &, + std::vector &, + std::vector &, + const bool); private: - int *nints; - void generate_fclists(int); + bool is_ascending(const int, const int *); }; } +// Define a hash function for FcProperty class +// Use boost::hash_combine +namespace std +{ + template <> + struct hash + { + std::size_t operator ()(ALM_NS::FcProperty const &obj) const + { + hash hasher; + size_t seed = 0; + for (auto i : obj.elems) { + seed ^= hasher(i) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + } + return seed; + } + }; +} diff --git a/alm/fitting.cpp b/alm/fitting.cpp index 61460e17..28fae7fa 100644 --- a/alm/fitting.cpp +++ b/alm/fitting.cpp @@ -67,7 +67,7 @@ void Fitting::fitmain() { int i; int nat = system->nat; - int natmin = symmetry->natmin; + int natmin = symmetry->nat_prim; int ntran = symmetry->ntran; int ndata = system->ndata; @@ -113,7 +113,7 @@ void Fitting::fitmain() N = 0; for (i = 0; i < maxorder; ++i) { - N += fcs->ndup[i].size(); + N += fcs->nequiv[i].size(); } std::cout << " Total Number of Parameters : " << N << std::endl << std::endl; @@ -124,7 +124,6 @@ void Fitting::fitmain() if (constraint->constraint_algebraic) { - N_new = 0; for (i = 0; i < maxorder; ++i) { N_new += constraint->index_bimap[i].size(); @@ -133,10 +132,10 @@ void Fitting::fitmain() << N_new << std::endl << std::endl; // memory->allocate(amat, M, N_new); - memory->allocate(amat_1D, N_new * M); + unsigned long NM = static_cast(N_new) * static_cast(M); + memory->allocate(amat_1D, NM); memory->allocate(fsum, M); memory->allocate(fsum_orig, M); - calc_matrix_elements_algebraic_constraint(M, N, N_new, nat, natmin, ndata_used, nmulti, maxorder, u, f, amat_1D, fsum, fsum_orig); @@ -379,8 +378,8 @@ void Fitting::data_multiplier(const int nat, f_rot[k] = f_tmp[3 * nat * i + 3 * j + k]; } - rotvec(u_rot, u_rot, symmetry->symrel[isym]); - rotvec(f_rot, f_rot, symmetry->symrel[isym]); + rotvec(u_rot, u_rot, symmetry->SymmData[isym].rotation_cart); + rotvec(f_rot, f_rot, symmetry->SymmData[isym].rotation_cart); for (k = 0; k < 3; ++k) { u[nmulti * idata + isym][3 * n_mapped + k] = u_rot[k]; @@ -710,7 +709,7 @@ void Fitting::fit_algebraic_constraints(int N, param_out[constraint->const_relate[i][j].p_index_target + ishift] = -tmp; } - ishift += fcs->ndup[i].size(); + ishift += fcs->nequiv[i].size(); iparam += constraint->index_bimap[i].size(); } @@ -780,7 +779,7 @@ void Fitting::fit_bootstrap(int N, ofs_fcs_boot << "# Relative Error(%), FCs ..."; for (i = 0; i < interaction->maxorder; ++i) { - ofs_fcs_boot << std::setw(10) << fcs->ndup[i].size(); + ofs_fcs_boot << std::setw(10) << fcs->nequiv[i].size(); } ofs_fcs_boot << std::endl; @@ -916,7 +915,7 @@ void Fitting::fit_consecutively(int N, ofs_fcs_seq << "# Relative Error(%), FCS..."; for (i = 0; i < interaction->maxorder; ++i) { - ofs_fcs_seq << std::setw(10) << fcs->ndup[i].size(); + ofs_fcs_seq << std::setw(10) << fcs->nequiv[i].size(); } ofs_fcs_seq << std::endl; @@ -1053,22 +1052,23 @@ void Fitting::calc_matrix_elements(const int M, mm = 0; - for (std::vector::iterator iter = fcs->ndup[order].begin(); - iter != fcs->ndup[order].end(); ++iter) { + for (std::vector::iterator iter = fcs->nequiv[order].begin(); + iter != fcs->nequiv[order].end(); ++iter) { for (i = 0; i < *iter; ++i) { - ind[0] = fcs->fc_set[order][mm].elems[0]; - k = idata + inprim_index(fcs->fc_set[order][mm].elems[0]); + ind[0] = fcs->fc_table[order][mm].elems[0]; + k = idata + inprim_index(fcs->fc_table[order][mm].elems[0]); amat_tmp = 1.0; for (j = 1; j < order + 2; ++j) { - ind[j] = fcs->fc_set[order][mm].elems[j]; - amat_tmp *= u[irow][fcs->fc_set[order][mm].elems[j]]; + ind[j] = fcs->fc_table[order][mm].elems[j]; + amat_tmp *= u[irow][fcs->fc_table[order][mm].elems[j]]; } - amat[k][iparam] -= gamma(order + 2, ind) * fcs->fc_set[order][mm].coef * amat_tmp; + amat[k][iparam] -= gamma(order + 2, ind) * fcs->fc_table[order][mm].sign * amat_tmp; ++mm; } ++iparam; } } + } memory->deallocate(ind); @@ -1093,15 +1093,14 @@ void Fitting::calc_matrix_elements_algebraic_constraint(const int M, double *bvec, double *bvec_orig) { - int i, j; - int irow; - int ncycle; + long i, j; + long irow; + long ncycle; std::cout << " Calculation of matrix elements for direct fitting started ... "; ncycle = ndata_fit * nmulti; - int natmin3 = 3 * natmin; - + long natmin3 = 3 * static_cast(natmin); #ifdef _OPENMP #pragma omp parallel for private(j) @@ -1119,10 +1118,10 @@ void Fitting::calc_matrix_elements_algebraic_constraint(const int M, #endif { int *ind; - int mm, order, iat, k; - int im, idata, iparam; - int ishift; - int iold, inew; + long mm, order, iat, k; + long im, idata, iparam; + long ishift; + long iold, inew; double amat_tmp; double **amat_orig; double **amat_mod; @@ -1164,18 +1163,18 @@ void Fitting::calc_matrix_elements_algebraic_constraint(const int M, mm = 0; - for (std::vector::iterator iter = fcs->ndup[order].begin(); - iter != fcs->ndup[order].end(); ++iter) { + for (std::vector::iterator iter = fcs->nequiv[order].begin(); + iter != fcs->nequiv[order].end(); ++iter) { for (i = 0; i < *iter; ++i) { - ind[0] = fcs->fc_set[order][mm].elems[0]; + ind[0] = fcs->fc_table[order][mm].elems[0]; k = inprim_index(ind[0]); amat_tmp = 1.0; for (j = 1; j < order + 2; ++j) { - ind[j] = fcs->fc_set[order][mm].elems[j]; - amat_tmp *= u[irow][fcs->fc_set[order][mm].elems[j]]; + ind[j] = fcs->fc_table[order][mm].elems[j]; + amat_tmp *= u[irow][fcs->fc_table[order][mm].elems[j]]; } - amat_orig[k][iparam] -= gamma(order + 2, ind) * fcs->fc_set[order][mm].coef * amat_tmp; + amat_orig[k][iparam] -= gamma(order + 2, ind) * fcs->fc_table[order][mm].sign * amat_tmp; ++mm; } ++iparam; @@ -1220,7 +1219,7 @@ void Fitting::calc_matrix_elements_algebraic_constraint(const int M, } } - ishift += fcs->ndup[order].size(); + ishift += fcs->nequiv[order].size(); iparam += constraint->index_bimap[order].size(); } @@ -1231,15 +1230,12 @@ void Fitting::calc_matrix_elements_algebraic_constraint(const int M, amat[natmin3 * ncycle * j + i + idata] = amat_mod[i][j]; } } - } memory->deallocate(ind); memory->deallocate(amat_orig); memory->deallocate(amat_mod); } - - std::cout << "done!" << std::endl << std::endl; } @@ -1249,7 +1245,7 @@ int Fitting::inprim_index(const int n) int atmn = n / 3; int crdn = n % 3; - for (int i = 0; i < symmetry->natmin; ++i) { + for (int i = 0; i < symmetry->nat_prim; ++i) { if (symmetry->map_p2s[i][0] == atmn) { in = 3 * i + crdn; break; diff --git a/alm/input.cpp b/alm/input.cpp index ab0a3eaf..3b4853b2 100644 --- a/alm/input.cpp +++ b/alm/input.cpp @@ -107,10 +107,11 @@ void Input::parse_general_vars() std::string *kdname; double **magmom, magmag; double tolerance; + double tolerance_constraint; std::vector kdname_v, periodic_v, magmom_v, str_split; std::string str_allowed_list = "PREFIX MODE NAT NKD NSYM KD PERIODIC PRINTSYM TOLERANCE DBASIS TRIMEVEN\ - MAGMOM NONCOLLINEAR TREVSYM HESSIAN"; + MAGMOM NONCOLLINEAR TREVSYM HESSIAN TOL_CONST"; std::string str_no_defaults = "PREFIX MODE NAT NKD KD"; std::vector no_defaults; std::map general_var_dict; @@ -197,6 +198,12 @@ void Input::parse_general_vars() assign_val(tolerance, "TOLERANCE", general_var_dict); } + if (general_var_dict["TOL_CONST"].empty()) { + tolerance_constraint = eps6; + } else { + assign_val(tolerance_constraint, "TOL_CONST", general_var_dict); + } + // Convert MAGMOM input to array memory->allocate(magmom, nat, 3); lspin = false; @@ -342,6 +349,7 @@ void Input::parse_general_vars() system->noncollinear = noncollinear; symmetry->trev_sym_mag = trevsym; writes->print_hessian = print_hessian; + constraint->tolerance_constraint = tolerance_constraint; if (mode == "suggest") { displace->disp_basis = str_disp_basis; @@ -449,6 +457,8 @@ void Input::parse_cell_parameter() system->lavec[i][j] = a * lavec_tmp[i][j]; } } + line_vec.clear(); + line_split.clear(); } void Input::parse_interaction_vars() diff --git a/alm/interaction.cpp b/alm/interaction.cpp index b9137c26..ae429b08 100644 --- a/alm/interaction.cpp +++ b/alm/interaction.cpp @@ -30,6 +30,7 @@ using namespace ALM_NS; Interaction::Interaction(ALM *alm) : Pointers(alm) { + rcs = nullptr; } Interaction::~Interaction() @@ -43,6 +44,9 @@ Interaction::~Interaction() memory->deallocate(interaction_pair); memory->deallocate(mindist_cluster); memory->deallocate(distall); + if (rcs) { + memory->deallocate(rcs); + } } void Interaction::init() @@ -79,8 +83,8 @@ void Interaction::init() memory->allocate(exist_image, nneib); memory->allocate(distall, nat, nat); memory->allocate(mindist_pairs, nat, nat); - memory->allocate(interaction_pair, maxorder, symmetry->natmin); - memory->allocate(mindist_cluster, maxorder, symmetry->natmin); + memory->allocate(interaction_pair, maxorder, symmetry->nat_prim); + memory->allocate(mindist_cluster, maxorder, symmetry->nat_prim); memory->allocate(pairs, maxorder); generate_coordinate_of_periodic_images(nat, system->xcoord, @@ -105,7 +109,7 @@ void Interaction::generate_pairs(std::set *pair_out, int i, j; int iat; int order; - int natmin = symmetry->natmin; + int natmin = symmetry->nat_prim; int nat = system->nat; int *pair_tmp; @@ -193,6 +197,8 @@ void Interaction::generate_coordinate_of_periodic_images(const unsigned int nat, for (ja = -1; ja <= 1; ++ja) { for (ka = -1; ka <= 1; ++ka) { + if (ia == 0 && ja == 0 && ka == 0) continue; + ++icell; // When periodic flag is zero along an axis, @@ -291,9 +297,9 @@ void Interaction::print_neighborlist(std::vector **mindist) double dist_tmp; std::vector *neighborlist; - memory->allocate(neighborlist, symmetry->natmin); + memory->allocate(neighborlist, symmetry->nat_prim); - for (i = 0; i < symmetry->natmin; ++i) { + for (i = 0; i < symmetry->nat_prim; ++i) { neighborlist[i].clear(); iat = symmetry->map_p2s[i][0]; @@ -312,7 +318,7 @@ void Interaction::print_neighborlist(std::vector **mindist) int nthnearest; std::vector atomlist; - for (i = 0; i < symmetry->natmin; ++i) { + for (i = 0; i < symmetry->nat_prim; ++i) { nthnearest = 0; atomlist.clear(); @@ -400,7 +406,7 @@ void Interaction::search_interactions(std::vector **interaction_list_out, // int i; int order; - int natmin = symmetry->natmin; + int natmin = symmetry->nat_prim; int iat, jat; int nat = system->nat; int ikd, jkd; @@ -706,7 +712,7 @@ void Interaction::calc_mindist_clusters(std::vector **interaction_pair_in, // Calculate the complete set of interaction clusters for each order. // - int natmin = symmetry->natmin; + int natmin = symmetry->nat_prim; int i, j, k; int iat, jat; int order; @@ -737,7 +743,7 @@ void Interaction::calc_mindist_clusters(std::vector **interaction_pair_in, std::vector> data_vec; for (order = 0; order < maxorder; ++order) { - for (i = 0; i < symmetry->natmin; ++i) { + for (i = 0; i < symmetry->nat_prim; ++i) { mindist_cluster_out[order][i].clear(); @@ -973,7 +979,7 @@ void Interaction::calc_mindist_clusters2(std::vector **interaction_pair_in, { std::vector distance_list; - int natmin = symmetry->natmin; + int natmin = symmetry->nat_prim; int i, j, k; int iat, jat; int ikd, jkd; @@ -998,7 +1004,7 @@ void Interaction::calc_mindist_clusters2(std::vector **interaction_pair_in, bool isok; for (order = 0; order < maxorder; ++order) { - for (i = 0; i < symmetry->natmin; ++i) { + for (i = 0; i < symmetry->nat_prim; ++i) { mindist_cluster_out[order][i].clear(); iat = symmetry->map_p2s[i][0]; diff --git a/alm/interaction.h b/alm/interaction.h index a4a3036e..c62cef8a 100644 --- a/alm/interaction.h +++ b/alm/interaction.h @@ -177,13 +177,13 @@ namespace ALM_NS { public: std::vector atom; - std::vector > cell; + std::vector> cell; double distmax; MinimumDistanceCluster(); MinimumDistanceCluster(const std::vector atom_in, - const std::vector > cell_in, + const std::vector> cell_in, const double dist_in) { for (int i = 0; i < atom_in.size(); ++i) { @@ -196,7 +196,7 @@ namespace ALM_NS } MinimumDistanceCluster(const std::vector atom_in, - const std::vector > cell_in) + const std::vector> cell_in) { for (int i = 0; i < atom_in.size(); ++i) { atom.push_back(atom_in[i]); @@ -279,11 +279,10 @@ namespace ALM_NS std::vector **, int *, std::set **); - void cell_combination(std::vector >, + void cell_combination(std::vector>, int, std::vector, - std::vector > &); + std::vector> &); void generate_pairs(std::set *, std::set **); }; } - diff --git a/alm/memory.h b/alm/memory.h index 4da3cdc7..b03fd92e 100644 --- a/alm/memory.h +++ b/alm/memory.h @@ -25,8 +25,8 @@ namespace ALM_NS // allocator - template - T* allocate(T *&arr, const unsigned int n1) + template + T* allocate(T *&arr, const A n1) { try { arr = new T [n1]; @@ -144,9 +144,11 @@ namespace ALM_NS // memsize calculator - unsigned long memsize_in_MB(const int size_of_one, const unsigned int n1) + + template + unsigned long memsize_in_MB(const int size_of_one, const A n1) { - unsigned long n = n1 * size_of_one; + unsigned long n = static_cast(n1) * size_of_one; return n / 1000000; } @@ -169,4 +171,3 @@ namespace ALM_NS } }; } - diff --git a/alm/patterndisp.cpp b/alm/patterndisp.cpp index 4cebd8b0..5a832028 100644 --- a/alm/patterndisp.cpp +++ b/alm/patterndisp.cpp @@ -40,6 +40,7 @@ void Displace::gen_displacement_pattern() { int i, j, m, order; int maxorder = interaction->maxorder; + std::string preferred_basis; std::vector group_tmp; std::vector *constsym; @@ -47,17 +48,50 @@ void Displace::gen_displacement_pattern() std::set *include_set; std::set *dispset; + std::vector *nequiv; + std::vector *fc_table, *fc_zeros; + std::vector *const_fix_tmp; std::vector *const_relate_tmp; boost::bimap *index_bimap_tmp; + std::cout << " DISPLACEMENT PATTERN" << std::endl; + std::cout << " ====================" << std::endl << std::endl; + + // Decide preferred basis (Cartesian or Lattice) + int ncompat_cart = 0; + int ncompat_latt = 0; + for (auto it = symmetry->SymmData.begin(); it != symmetry->SymmData.end(); ++it) { + if ((*it).compatible_with_cartesian) ++ncompat_cart; + if ((*it).compatible_with_lattice) ++ncompat_latt; + } + if (ncompat_cart >= ncompat_latt) { + preferred_basis = "Cartesian"; + } else { + preferred_basis = "Lattice"; + } + + memory->allocate(fc_table, maxorder); + memory->allocate(fc_zeros, maxorder); memory->allocate(constsym, maxorder); + memory->allocate(nequiv, maxorder); memory->allocate(const_fix_tmp, maxorder); memory->allocate(const_relate_tmp, maxorder); memory->allocate(index_bimap_tmp, maxorder); - constraint->constraint_from_symmetry(constsym); - constraint->get_mapping_constraint(maxorder, constsym, const_fix_tmp, + for (order = 0; order < maxorder; ++order) { + fcs->generate_force_constant_table(order, interaction->pairs[order], + symmetry->SymmData, preferred_basis, + fc_table[order], nequiv[order], + fc_zeros[order], false); + + constraint->get_symmetry_constraint(order, interaction->pairs[order], + symmetry->SymmData, preferred_basis, + fc_table[order], nequiv[order], + constsym[order]); + } + + constraint->get_mapping_constraint(maxorder, nequiv, constsym, const_fix_tmp, const_relate_tmp, index_bimap_tmp, true); for (order = 0; order < maxorder; ++order) { @@ -70,6 +104,7 @@ void Displace::gen_displacement_pattern() memory->deallocate(constsym); memory->deallocate(const_fix_tmp); memory->deallocate(const_relate_tmp); + memory->deallocate(fc_zeros); memory->allocate(include_set, maxorder); @@ -85,12 +120,12 @@ void Displace::gen_displacement_pattern() memory->deallocate(index_bimap_tmp); - std::cout << " Generating displacement patterns in "; - if (disp_basis[0] == 'C') { - std::cout << "Cartesian coordinate..."; - } else { - std::cout << "fractional coordinate..."; - } + std::cout << " Generating displacement patterns in "; + // if (disp_basis[0] == 'C') { + std::cout << "Cartesian coordinate... "; + // } else { + // std::cout << "fractional coordinate..."; + // } memory->allocate(dispset, maxorder); @@ -98,7 +133,7 @@ void Displace::gen_displacement_pattern() m = 0; - for (i = 0; i < fcs->ndup[order].size(); ++i) { + for (i = 0; i < nequiv[order].size(); ++i) { if (include_set[order].find(i) != include_set[order].end()) { @@ -108,7 +143,7 @@ void Displace::gen_displacement_pattern() // Here, duplicate entries will be removed. // For example, (iij) will be reduced to (ij). for (j = 0; j < order + 1; ++j) { - group_tmp.push_back(fcs->fc_set[order][m].elems[j]); + group_tmp.push_back(fc_table[order][m].elems[j]); } group_tmp.erase(std::unique(group_tmp.begin(), group_tmp.end()), group_tmp.end()); @@ -118,23 +153,35 @@ void Displace::gen_displacement_pattern() } - m += fcs->ndup[order][i]; + m += nequiv[order][i]; } } memory->deallocate(include_set); + memory->deallocate(nequiv); + memory->deallocate(fc_table); memory->allocate(pattern_all, maxorder); - generate_pattern_all(maxorder, pattern_all, dispset); + generate_pattern_all(maxorder, pattern_all, + dispset, preferred_basis); memory->deallocate(dispset); std::cout << " done!" << std::endl; + std::cout << std::endl; + + for (order = 0; order < maxorder; ++order) { + std::cout << " Number of disp. patterns for " << std::setw(9) + << interaction->str_order[order] << " : " + << pattern_all[order].size() << std::endl; + } + std::cout << std::endl; } void Displace::generate_pattern_all(const int N, std::vector *pattern, - std::set *dispset_in) + std::set *dispset_in, + const std::string preferred_basis) { int i, j; int order; @@ -158,8 +205,7 @@ void Displace::generate_pattern_all(const int N, pattern[order].clear(); - for (std::set::iterator it = dispset_in[order].begin(); - it != dispset_in[order].end(); ++it) { + for (auto it = dispset_in[order].begin(); it != dispset_in[order].end(); ++it) { atoms.clear(); directions.clear(); @@ -184,7 +230,7 @@ void Displace::generate_pattern_all(const int N, if (trim_dispsign_for_evenfunc) { find_unique_sign_pairs(natom_disp, sign_prod[natom_disp - 1], - nums, sign_reduced); + nums, sign_reduced, preferred_basis); } else { sign_reduced.clear(); std::copy(sign_prod[natom_disp - 1].begin(), @@ -196,8 +242,7 @@ void Displace::generate_pattern_all(const int N, std::copy(directions.begin(), directions.end(), std::back_inserter(directions_copy)); - for (std::vector>::const_iterator it2 = sign_reduced.begin(); - it2 != sign_reduced.end(); ++it2) { + for (auto it2 = sign_reduced.cbegin(); it2 != sign_reduced.cend(); ++it2) { directions.clear(); for (i = 0; i < (*it2).size(); ++i) { @@ -207,11 +252,19 @@ void Displace::generate_pattern_all(const int N, disp_tmp[j] = directions_copy[3 * i + j] * sign_double; } - if (disp_basis[0] == 'F') { - rotvec(disp_tmp, disp_tmp, system->rlavec); - for (j = 0; j < 3; ++j) { - disp_tmp[j] /= 2.0 * pi; - } + // if (disp_basis[0] == 'F') { + // rotvec(disp_tmp, disp_tmp, system->rlavec); + // for (j = 0; j < 3; ++j) { + // disp_tmp[j] /= 2.0 * pi; + // } + // } + + if (preferred_basis == "Lattice") { + rotvec(disp_tmp, disp_tmp, system->lavec); + double norm = 0.0; + for (j = 0; j < 3; ++j) norm += disp_tmp[j] * disp_tmp[j]; + norm = std::sqrt(norm); + for (j = 0; j < 3; ++j) disp_tmp[j] /= norm; } for (j = 0; j < 3; ++j) { @@ -254,7 +307,8 @@ void Displace::generate_signvecs(const int N, void Displace::find_unique_sign_pairs(const int N, std::vector> sign_in, std::vector pair_in, - std::vector> &sign_out) + std::vector> &sign_out, + const std::string preferred_basis) { int isym, i, j, k; int mapped_atom; @@ -334,9 +388,21 @@ void Displace::find_unique_sign_pairs(const int N, for (j = 0; j < 3; ++j) { disp_sym[mapped_atom][j] = 0.0; - for (k = 0; k < 3; ++k) { - disp_sym[mapped_atom][j] - += symmetry->symrel[isym][j][k] * disp[list_disp_atom[i]][k]; + + if (preferred_basis == "Cartesian") { + for (k = 0; k < 3; ++k) { + disp_sym[mapped_atom][j] + += symmetry->SymmData[isym].rotation_cart[j][k] * disp[list_disp_atom[i]][k]; + } + } else if (preferred_basis == "Lattice") { + for (k = 0; k < 3; ++k) { + disp_sym[mapped_atom][j] + += static_cast(symmetry->SymmData[isym].rotation[j][k]) + * disp[list_disp_atom[i]][k]; + } + } else { + error->exit("find_unique_sign_pairs", + "Invalid basis. This cannot happen."); } disp_tmp = disp_sym[mapped_atom][j]; @@ -388,10 +454,23 @@ void Displace::find_unique_sign_pairs(const int N, for (j = 0; j < 3; ++j) { disp_sym[mapped_atom][j] = 0.0; - for (k = 0; k < 3; ++k) { - disp_sym[mapped_atom][j] += symmetry->symrel[symnum_vec[isym]][j][k] - * disp[list_disp_atom[i]][k]; + if (preferred_basis == "Cartesian") { + for (k = 0; k < 3; ++k) { + disp_sym[mapped_atom][j] + += symmetry->SymmData[symnum_vec[isym]].rotation_cart[j][k] + * disp[list_disp_atom[i]][k]; + } + } else if (preferred_basis == "Lattice") { + for (k = 0; k < 3; ++k) { + disp_sym[mapped_atom][j] + += static_cast(symmetry->SymmData[symnum_vec[isym]].rotation[j][k]) + * disp[list_disp_atom[i]][k]; + } + } else { + error->exit("find_unique_sign_pairs", + "Invalid basis. This cannot happen."); } + disp_tmp = disp_sym[mapped_atom][j]; if (std::abs(disp_tmp) > eps) { diff --git a/alm/patterndisp.h b/alm/patterndisp.h index 768dd3e8..387acaa7 100644 --- a/alm/patterndisp.h +++ b/alm/patterndisp.h @@ -26,7 +26,7 @@ namespace ALM_NS DispAtomSet(std::vector vec) { - for (std::vector::iterator it = vec.begin(); it != vec.end(); ++it) { + for (auto it = vec.begin(); it != vec.end(); ++it) { atomset.push_back((*it)); } } @@ -56,7 +56,7 @@ namespace ALM_NS DispDirectionHarmonic(int n, std::vector list_in) { atom = n; - for (std::vector::iterator it = list_in.begin(); it != list_in.end(); ++it) { + for (auto it = list_in.begin(); it != list_in.end(); ++it) { directionlist.push_back(*it); } } @@ -134,15 +134,17 @@ namespace ALM_NS std::vector disp_harm, disp_harm_best; void generate_pattern_all(const int, std::vector *, - std::set *); + std::set *, + const std::string); void generate_signvecs(const int, - std::vector > &, + std::vector> &, std::vector); void find_unique_sign_pairs(const int, - std::vector >, + std::vector>, std::vector, - std::vector > &); + std::vector> &, + const std::string); }; } diff --git a/alm/pointers.h b/alm/pointers.h index f1bec298..77e7af4b 100644 --- a/alm/pointers.h +++ b/alm/pointers.h @@ -31,9 +31,13 @@ namespace ALM_NS displace(ptr->displace), writes(ptr->writes), error(ptr->error), - timer(ptr->timer) {} + timer(ptr->timer) + { + } - virtual ~Pointers() {} + virtual ~Pointers() + { + } protected: ALM *alm; @@ -52,4 +56,3 @@ namespace ALM_NS Timer *&timer; }; } - diff --git a/alm/symmetry.cpp b/alm/symmetry.cpp index d9b13875..a78a7c30 100644 --- a/alm/symmetry.cpp +++ b/alm/symmetry.cpp @@ -31,19 +31,15 @@ using namespace ALM_NS; Symmetry::Symmetry(ALM *alm) : Pointers(alm) { - file_sym = "SYMM_INFO"; } Symmetry::~Symmetry() { - memory->deallocate(symrel); - memory->deallocate(symrel_int); - memory->deallocate(tnons); memory->deallocate(map_sym); memory->deallocate(map_p2s); memory->deallocate(map_s2p); memory->deallocate(symnum_tran); - memory->deallocate(sym_available); + SymmData.clear(); } void Symmetry::init() @@ -54,49 +50,48 @@ void Symmetry::init() std::cout << " SYMMETRY" << std::endl; std::cout << " ========" << std::endl << std::endl; - setup_symmetry_operation(nat, nsym, system->lavec, system->rlavec, + setup_symmetry_operation(nat, nsym, + system->lavec, system->rlavec, system->xcoord, system->kd); - memory->allocate(tnons, nsym, 3); - memory->allocate(symrel_int, nsym, 3, 3); - - int isym = 0; - for (std::vector::iterator iter = SymmList.begin(); - iter != SymmList.end(); ++iter) { - for (i = 0; i < 3; ++i) { - for (j = 0; j < 3; ++j) { - symrel_int[isym][i][j] = (*iter).rot[i][j]; - } - } - for (i = 0; i < 3; ++i) { - tnons[isym][i] = (*iter).tran[i]; - } - ++isym; - } - - std::cout << " Number of symmetry operations = " << nsym << std::endl; - memory->allocate(symrel, nsym, 3, 3); - symop_in_cart(system->lavec, system->rlavec); - - memory->allocate(sym_available, nsym); - int nsym_fc; - symop_availability_check(symrel, sym_available, nsym, nsym_fc); - - if (nsym_fc == nsym) { - std::cout << " All symmetry operations will be used to" << std::endl; - std::cout << " reduce the number of force constants." << std::endl; - } else { - std::cout << " " << nsym_fc << " symmetry operations out of " - << nsym << " will be used to reduce the number of parameters." << std::endl; - std::cout << " Other " << nsym - nsym_fc - << " symmetry operations will be imposed as constraints." << std::endl; - } - std::cout << std::endl; + std::cout << " Number of symmetry operations = " << SymmData.size() << std::endl; + + // int nsym_fc; + + // if (nsym_fc == nsym) { + // std::cout << " All symmetry operations will be used to" << std::endl; + // std::cout << " reduce the number of force constants." << std::endl; + // } else { + // std::cout << " " << nsym_fc << " symmetry operations out of " + // << nsym << " will be used to reduce the number of parameters." << std::endl; + // std::cout << " Other " << nsym - nsym_fc + // << " symmetry operations will be imposed as constraints." << std::endl; + // } + // std::cout << std::endl; + + // int counter = 0; + // for (auto it = SymmData.begin(); it != SymmData.end(); ++it) { + // std::cout << "Symm. No. : " << std::setw(4) << counter + 1; + // std::cout << "( " << (*it).compatible_with_lattice << " " << (*it).compatible_with_cartesian << ")" << std::endl; + // for (i = 0; i < 3; ++i) { + // for (j = 0; j < 3; ++j) { + // std::cout << std::setw(15) << (*it).rotation[i][j]; + // } + // std::cout << " "; + // for (j = 0; j < 3; ++j) { + // std::cout << std::setw(15) << (*it).rotation_cart[i][j]; + // } + // + // std::cout << std::endl; + // } + // std::cout << std::endl; + // ++counter; + // } pure_translations(); memory->allocate(map_sym, nat, nsym); - memory->allocate(map_p2s, natmin, ntran); + memory->allocate(map_p2s, nat_prim, ntran); memory->allocate(map_s2p, nat); genmaps(nat, system->xcoord, map_sym, map_p2s, map_s2p); @@ -107,7 +102,7 @@ void Symmetry::init() for (int i = 0; i < ntran; ++i) { std::cout << std::setw(6) << i + 1 << " | "; - for (int j = 0; j < natmin; ++j) { + for (int j = 0; j < nat_prim; ++j) { std::cout << std::setw(5) << map_p2s[j][i] + 1; if ((j + 1) % 5 == 0) { std::cout << std::endl << " | "; @@ -131,7 +126,7 @@ void Symmetry::setup_symmetry_operation(int nat, { int i, j; - SymmList.clear(); + SymmData.clear(); if (nsym == 0) { @@ -141,11 +136,11 @@ void Symmetry::setup_symmetry_operation(int nat, std::cout << " Please be patient. " << std::endl; std::cout << " This can take a while for a large supercell." << std::endl << std::endl; - findsym(nat, aa, x, SymmList); - // The order in SymmList changes for each run because it was generated + findsym(nat, aa, x, SymmData); + // The order in SymmData changes for each run because it was generated // with OpenMP. Therefore, we sort the list here to have the same result. - std::sort(SymmList.begin() + 1, SymmList.end()); - nsym = SymmList.size(); + std::sort(SymmData.begin() + 1, SymmData.end()); + nsym = SymmData.size(); if (is_printsymmetry) { std::ofstream ofs_sym; @@ -154,16 +149,15 @@ void Symmetry::setup_symmetry_operation(int nat, ofs_sym.open(file_sym.c_str(), std::ios::out); ofs_sym << nsym << std::endl; - for (std::vector::iterator p = SymmList.begin(); - p != SymmList.end(); ++p) { + for (auto p = SymmData.begin(); p != SymmData.end(); ++p) { for (i = 0; i < 3; ++i) { for (j = 0; j < 3; ++j) { - ofs_sym << std::setw(4) << (*p).rot[i][j]; + ofs_sym << std::setw(4) << (*p).rotation[i][j]; } } ofs_sym << " "; for (i = 0; i < 3; ++i) { - ofs_sym << std::setprecision(15) << std::setw(20) << (*p).tran[i]; + ofs_sym << std::setprecision(15) << std::setw(21) << (*p).tran[i]; } ofs_sym << std::endl; } @@ -175,30 +169,41 @@ void Symmetry::setup_symmetry_operation(int nat, // Identity operation only ! - std::cout << " NSYM = 1 : Only the identity matrix will be considered." << std::endl << std::endl; + std::cout << " NSYM = 1 : Only the identity matrix will be considered." + << std::endl << std::endl; int rot_tmp[3][3]; + double rot_cart_tmp[3][3]; double tran_tmp[3]; for (i = 0; i < 3; ++i) { for (j = 0; j < 3; ++j) { if (i == j) { rot_tmp[i][j] = 1; + rot_cart_tmp[i][j] = 1.0; } else { rot_tmp[i][j] = 0; + rot_cart_tmp[i][j] = 0.0; } } tran_tmp[i] = 0.0; } - SymmList.push_back(SymmetryOperation(rot_tmp, tran_tmp)); + SymmData.push_back(SymmetryOperation(rot_tmp, + tran_tmp, + rot_cart_tmp, + true, + true, + true)); } else { - std::cout << " NSYM > 1 : Symmetry operations will be read from SYMM_INFO file" << std::endl << std::endl; + std::cout << " NSYM > 1 : Symmetry operations will be read from SYMM_INFO file" + << std::endl << std::endl; int nsym2; int rot_tmp[3][3]; + double rot_cart_tmp[3][3]; double tran_tmp[3]; std::ifstream ifs_sym; @@ -216,13 +221,19 @@ void Symmetry::setup_symmetry_operation(int nat, >> rot_tmp[2][0] >> rot_tmp[2][1] >> rot_tmp[2][2] >> tran_tmp[0] >> tran_tmp[1] >> tran_tmp[2]; - SymmList.push_back(SymmetryOperation(rot_tmp, tran_tmp)); + symop_in_cart(rot_cart_tmp, rot_tmp, system->lavec, system->rlavec); + SymmData.push_back(SymmetryOperation(rot_tmp, + tran_tmp, + rot_cart_tmp, + is_compatible(rot_tmp), + is_compatible(rot_cart_tmp), + is_translation(rot_tmp))); } ifs_sym.close(); } #ifdef _DEBUG - print_symmetrized_coordinate(x); + // print_symmetrized_coordinate(x); #endif } @@ -402,18 +413,24 @@ void Symmetry::find_crystal_symmetry(int nat, for (j = 0; j < 3; ++j) { if (i == j) { rot_int[i][j] = 1; + rot_cart[i][j] = 1.0; } else { rot_int[i][j] = 0; + rot_cart[i][j] = 0.0; } } tran[i] = 0.0; } - CrystalSymmList.push_back(SymmetryOperation(rot_int, tran)); + CrystalSymmList.push_back(SymmetryOperation(rot_int, + tran, + rot_cart, + is_compatible(rot_int), + is_compatible(rot_cart), + is_translation(rot_int))); - for (std::vector::iterator it_latsym = LatticeSymmList.begin(); - it_latsym != LatticeSymmList.end(); ++it_latsym) { + for (auto it_latsym = LatticeSymmList.begin(); it_latsym != LatticeSymmList.end(); ++it_latsym) { iat = atomclass[0][0]; @@ -484,12 +501,7 @@ void Symmetry::find_crystal_symmetry(int nat, } } - if (isok && system->lspin && system->noncollinear) { - for (i = 0; i < 3; ++i) { - mag[i] = system->magmom[jat][i]; - mag_rot[i] = system->magmom[iat][i]; - } - + if (isok) { matmul3(rot_tmp, rot, system->rlavec); matmul3(rot_cart, system->lavec, rot_tmp); @@ -498,36 +510,51 @@ void Symmetry::find_crystal_symmetry(int nat, rot_cart[i][j] /= (2.0 * pi); } } - rotvec(mag_rot, mag_rot, rot_cart); - // In the case of improper rotation, the factor -1 should be multiplied - // because the inversion operation doesn't flip the spin. - if (!is_proper(rot_cart)) { + if (system->lspin && system->noncollinear) { for (i = 0; i < 3; ++i) { - mag_rot[i] = -mag_rot[i]; + mag[i] = system->magmom[jat][i]; + mag_rot[i] = system->magmom[iat][i]; + } + + + rotvec(mag_rot, mag_rot, rot_cart); + + // In the case of improper rotation, the factor -1 should be multiplied + // because the inversion operation doesn't flip the spin. + if (!is_proper(rot_cart)) { + for (i = 0; i < 3; ++i) { + mag_rot[i] = -mag_rot[i]; + } } - } - mag_sym1 = (std::pow(mag[0] - mag_rot[0], 2.0) - + std::pow(mag[1] - mag_rot[1], 2.0) - + std::pow(mag[2] - mag_rot[2], 2.0)) < eps6; + mag_sym1 = (std::pow(mag[0] - mag_rot[0], 2.0) + + std::pow(mag[1] - mag_rot[1], 2.0) + + std::pow(mag[2] - mag_rot[2], 2.0)) < eps6; - mag_sym2 = (std::pow(mag[0] + mag_rot[0], 2.0) - + std::pow(mag[1] + mag_rot[1], 2.0) - + std::pow(mag[2] + mag_rot[2], 2.0)) < eps6; + mag_sym2 = (std::pow(mag[0] + mag_rot[0], 2.0) + + std::pow(mag[1] + mag_rot[1], 2.0) + + std::pow(mag[2] + mag_rot[2], 2.0)) < eps6; - if (!mag_sym1 && !mag_sym2) { - isok = false; - } else if (!mag_sym1 && mag_sym2 && !trev_sym_mag) { - isok = false; + if (!mag_sym1 && !mag_sym2) { + isok = false; + } else if (!mag_sym1 && mag_sym2 && !trev_sym_mag) { + isok = false; + } } } + if (isok) { #ifdef _OPENMP #pragma omp critical #endif - CrystalSymmList.push_back(SymmetryOperation((*it_latsym).mat, tran)); + CrystalSymmList.push_back(SymmetryOperation((*it_latsym).mat, + tran, + rot_cart, + is_compatible((*it_latsym).mat), + is_compatible(rot_cart), + is_translation((*it_latsym).mat))); } } @@ -535,54 +562,42 @@ void Symmetry::find_crystal_symmetry(int nat, } -void Symmetry::symop_in_cart(double lavec[3][3], double rlavec[3][3]) +void Symmetry::symop_in_cart(double rot_cart[3][3], + const int rot_lattice[3][3], + const double lavec[3][3], + const double rlavec[3][3]) { int i, j; - double sym_tmp[3][3], sym_crt[3][3]; + double sym_tmp[3][3]; double tmp[3][3]; - for (int isym = 0; isym < nsym; ++isym) { - - for (i = 0; i < 3; ++i) { - for (j = 0; j < 3; ++j) { - sym_tmp[i][j] = static_cast(symrel_int[isym][i][j]); - } - } - - matmul3(tmp, sym_tmp, rlavec); - matmul3(sym_crt, lavec, tmp); - - for (i = 0; i < 3; ++i) { - for (j = 0; j < 3; ++j) { - symrel[isym][i][j] = sym_crt[i][j] / (2.0 * pi); - } + for (i = 0; i < 3; ++i) { + for (j = 0; j < 3; ++j) { + sym_tmp[i][j] = static_cast(rot_lattice[i][j]); } } -#ifdef _DEBUG + matmul3(tmp, sym_tmp, rlavec); + matmul3(rot_cart, lavec, tmp); - std::cout << "Symmetry Operations in Cartesian Coordinate" << std::endl; - for (int isym = 0; isym < nsym; ++isym) { - for (i = 0; i < 3; ++i) { - for (j = 0; j < 3; ++j) { - std::cout << std::setw(8) << symrel[isym][i][j]; - } + for (i = 0; i < 3; ++i) { + for (j = 0; j < 3; ++j) { + rot_cart[i][j] = rot_cart[i][j] / (2.0 * pi); } - std::cout << std::endl; } -#endif } + void Symmetry::pure_translations() { int i; ntran = 0; for (i = 0; i < nsym; ++i) { - if (is_translation(symrel_int[i])) ++ntran; + if (SymmData[i].is_translation) ++ntran; } - natmin = system->nat / ntran; + nat_prim = system->nat / ntran; if (ntran > 1) { std::cout << " Given system is not primitive cell." << std::endl; @@ -591,11 +606,11 @@ void Symmetry::pure_translations() } else { std::cout << " Given system is a primitive cell." << std::endl; } - std::cout << " Primitive cell contains " << natmin << " atoms" << std::endl; + std::cout << " Primitive cell contains " << nat_prim << " atoms" << std::endl; if (system->nat % ntran) { error->exit("pure_translations", - "nat != natmin * ntran. Something is wrong in the structure."); + "nat != nat_prim * ntran. Something is wrong in the structure."); } memory->allocate(symnum_tran, ntran); @@ -603,7 +618,7 @@ void Symmetry::pure_translations() int isym = 0; for (i = 0; i < nsym; ++i) { - if (is_translation(symrel_int[i])) symnum_tran[isym++] = i; + if (SymmData[i].is_translation) symnum_tran[isym++] = i; } } @@ -634,7 +649,7 @@ void Symmetry::genmaps(int nat, for (i = 0; i < 3; ++i) { for (j = 0; j < 3; ++j) { - rot_double[i][j] = static_cast(symrel_int[isym][i][j]); + rot_double[i][j] = static_cast(SymmData[isym].rotation[i][j]); } } @@ -646,7 +661,7 @@ void Symmetry::genmaps(int nat, rotvec(xnew, x[iat], rot_double); - for (i = 0; i < 3; ++i) xnew[i] += tnons[isym][i]; + for (i = 0; i < 3; ++i) xnew[i] += SymmData[isym].tran[i]; for (jj = 0; jj < system->atomlist_class[itype].size(); ++jj) { @@ -691,7 +706,7 @@ void Symmetry::genmaps(int nat, memory->deallocate(is_checked); - for (iat = 0; iat < natmin; ++iat) { + for (iat = 0; iat < nat_prim; ++iat) { for (i = 0; i < ntran; ++i) { atomnum_translated = map_p2s[iat][i]; map_s2p[atomnum_translated].atom_num = iat; @@ -700,7 +715,7 @@ void Symmetry::genmaps(int nat, } } -bool Symmetry::is_translation(int **rot) +bool Symmetry::is_translation(const int rot[3][3]) { bool ret; @@ -712,32 +727,26 @@ bool Symmetry::is_translation(int **rot) return ret; } -void Symmetry::symop_availability_check(double ***rot, - bool *flag, - const int n, - int &nsym_fc) + +template +bool Symmetry::is_compatible(const T rot[3][3], + const double tolerance_zero) { - int i, j, k; + int i, j; int nfinite; + double rot_double[3][3]; - nsym_fc = 0; - - for (i = 0; i < nsym; ++i) { - - nfinite = 0; + nfinite = 0; + for (i = 0; i < 3; ++i) { for (j = 0; j < 3; ++j) { - for (k = 0; k < 3; ++k) { - if (std::abs(rot[i][j][k]) > 1.0e-5) ++nfinite; - } - } - - if (nfinite == 3) { - ++nsym_fc; - flag[i] = true; - } else { - flag[i] = false; + rot_double[i][j] = static_cast(rot[i][j]); + if (std::abs(rot_double[i][j]) > tolerance_zero) ++nfinite; } } + + if (nfinite == 3) return true; + + return false; } void Symmetry::print_symmetrized_coordinate(double **x) @@ -766,21 +775,21 @@ void Symmetry::print_symmetrized_coordinate(double **x) } } - for (std::vector::iterator it = SymmList.begin(); - it != SymmList.end(); ++it) { + for (std::vector::iterator it = SymmData.begin(); + it != SymmData.end(); ++it) { ++isym; std::cout << "Symmetry No. : " << std::setw(5) << isym << std::endl; - m11 = (*it).rot[0][0]; - m12 = (*it).rot[0][1]; - m13 = (*it).rot[0][2]; - m21 = (*it).rot[1][0]; - m22 = (*it).rot[1][1]; - m23 = (*it).rot[1][2]; - m31 = (*it).rot[2][0]; - m32 = (*it).rot[2][1]; - m33 = (*it).rot[2][2]; + m11 = (*it).rotation[0][0]; + m12 = (*it).rotation[0][1]; + m13 = (*it).rotation[0][2]; + m21 = (*it).rotation[1][0]; + m22 = (*it).rotation[1][1]; + m23 = (*it).rotation[1][2]; + m31 = (*it).rotation[2][0]; + m32 = (*it).rotation[2][1]; + m33 = (*it).rotation[2][2]; for (i = 0; i < 3; ++i) tran[i] = (*it).tran[i]; @@ -893,7 +902,7 @@ void Symmetry::print_symmetrized_coordinate(double **x) for (i = 0; i < nat; ++i) { for (j = 0; j < 3; ++j) { - x_avg[i][j] /= static_cast(SymmList.size()); + x_avg[i][j] /= static_cast(SymmData.size()); } } @@ -913,7 +922,7 @@ void Symmetry::print_symmetrized_coordinate(double **x) memory->deallocate(x_avg); } -bool Symmetry::is_proper(double rot[3][3]) +bool Symmetry::is_proper(const double rot[3][3]) { double det; bool ret; diff --git a/alm/symmetry.h b/alm/symmetry.h index 8da44287..33511447 100644 --- a/alm/symmetry.h +++ b/alm/symmetry.h @@ -1,7 +1,7 @@ /* symmetry.h - Copyright (c) 2014, 2015, 2016 Terumasa Tadano + Copyright (c) 2014--2017 Terumasa Tadano This file is distributed under the terms of the MIT license. Please see the file 'LICENCE.txt' in the root directory @@ -24,32 +24,44 @@ namespace ALM_NS class SymmetryOperation { public: - int rot[3][3]; - double tran[3]; + int rotation[3][3]; // in lattice basis + double tran[3]; // in Cartesian basis + double rotation_cart[3][3]; // in Cartesian basis + bool compatible_with_lattice; + bool compatible_with_cartesian; + bool is_translation; SymmetryOperation(); - // Declaration construction - - SymmetryOperation(const int rot_in[3][3], const double tran_in[3]) + SymmetryOperation(const int rot_in[3][3], + const double tran_in[3], + const double rot_cart_in[3][3], + const bool compatibility_lat, + const bool compatibility_cart, + const bool is_trans_in) { for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) { - rot[i][j] = rot_in[i][j]; + rotation[i][j] = rot_in[i][j]; + rotation_cart[i][j] = rot_cart_in[i][j]; } } for (int i = 0; i < 3; ++i) { tran[i] = tran_in[i]; } + compatible_with_lattice = compatibility_lat; + compatible_with_cartesian = compatibility_cart; + is_translation = is_trans_in; } + // Operator definition to sort bool operator<(const SymmetryOperation &a) const { std::vector v1, v2; for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) { - v1.push_back(static_cast(rot[i][j])); - v2.push_back(static_cast(a.rot[i][j])); + v1.push_back(static_cast(rotation[i][j])); + v2.push_back(static_cast(a.rotation[i][j])); } } for (int i = 0; i < 3; ++i) { @@ -94,14 +106,12 @@ namespace ALM_NS void init(); - unsigned int nsym, ntran, natmin; + unsigned int nsym, ntran, nat_prim; int is_printsymmetry; int multiply_data; int *symnum_tran; double tolerance; - double ***symrel; - double **tnons; int **map_sym; int **map_p2s; @@ -116,7 +126,8 @@ namespace ALM_NS Maps *map_s2p; int trev_sym_mag; - bool *sym_available; + + std::vector SymmData; private: @@ -130,24 +141,24 @@ namespace ALM_NS void findsym(int, double [3][3], double **, std::vector &); - bool is_translation(int **); - bool is_proper(double [3][3]); + bool is_translation(const int [3][3]); + bool is_proper(const double [3][3]); - void symop_in_cart(double [3][3], double [3][3]); + void symop_in_cart(double [3][3], const int [3][3], + const double [3][3], const double [3][3]); void pure_translations(); void print_symmetrized_coordinate(double **); - void symop_availability_check(double ***, bool *, const int, int &); + + template + bool is_compatible(const T [3][3], const double tolerance_zero = 1.0e-5); void find_lattice_symmetry(double [3][3], std::vector &); void find_crystal_symmetry(int, int, - std::vector *, double **x, + std::vector *, double **, std::vector, std::vector &); - std::string file_sym; - int ***symrel_int; - std::vector SymmList; + std::string file_sym = "SYMM_INFO"; }; } - diff --git a/alm/system.cpp b/alm/system.cpp index 70c85e05..3267c5f3 100644 --- a/alm/system.cpp +++ b/alm/system.cpp @@ -32,13 +32,34 @@ using namespace ALM_NS; System::System(ALM *alm): Pointers(alm) { + kdname = nullptr; + kd = nullptr; + xcoord = nullptr; + magmom = nullptr; + x_cartesian = nullptr; + atomlist_class = nullptr; } System::~System() { - memory->deallocate(x_cartesian); - memory->deallocate(atomlist_class); - memory->deallocate(magmom); + if (kdname) { + memory->deallocate(kdname); + } + if (kd) { + memory->deallocate(kd); + } + if (xcoord) { + memory->deallocate(xcoord); + } + if (magmom) { + memory->deallocate(magmom); + } + if (x_cartesian) { + memory->deallocate(x_cartesian); + } + if (atomlist_class) { + memory->deallocate(atomlist_class); + } } void System::init() @@ -251,7 +272,7 @@ void System::load_reference_system_xml(std::string file_reference_fcs, get_value_from_xml(pt, "Data.Symmetry.NumberOfTranslations")); natmin_ref = nat_ref / ntran_ref; - if (natmin_ref != symmetry->natmin) { + if (natmin_ref != symmetry->nat_prim) { error->exit("load_reference_system_xml", "The number of atoms in the primitive cell is not consistent."); } @@ -260,7 +281,7 @@ void System::load_reference_system_xml(std::string file_reference_fcs, nfcs_ref = boost::lexical_cast( get_value_from_xml(pt, "Data.ForceConstants.HarmonicUnique.NFC2")); - if (nfcs_ref != fcs->ndup[0].size()) { + if (nfcs_ref != fcs->nequiv[0].size()) { error->exit("load_reference_system_xml", "The number of harmonic force constants is not the same."); } @@ -269,7 +290,7 @@ void System::load_reference_system_xml(std::string file_reference_fcs, nfcs_ref = boost::lexical_cast( get_value_from_xml(pt, "Data.ForceConstants.CubicUnique.NFC3")); - if (nfcs_ref != fcs->ndup[1].size()) { + if (nfcs_ref != fcs->nequiv[1].size()) { error->exit("load_reference_system_xml", "The number of cubic force constants is not the same."); } @@ -316,13 +337,13 @@ void System::load_reference_system_xml(std::string file_reference_fcs, list_found.clear(); - for (std::vector::iterator p = fcs->fc_set[order_fcs].begin(); - p != fcs->fc_set[order_fcs].end(); ++p) { + for (std::vector::iterator p = fcs->fc_table[order_fcs].begin(); + p != fcs->fc_table[order_fcs].end(); ++p) { FcProperty list_tmp = *p; // Using copy constructor for (i = 0; i < nterms; ++i) { ind[i] = list_tmp.elems[i]; } - list_found.insert(FcProperty(nterms, list_tmp.coef, + list_found.insert(FcProperty(nterms, list_tmp.sign, ind, list_tmp.mother)); } @@ -368,7 +389,7 @@ void System::load_reference_system() bool is_found_system = false; int nparam_harmonic_ref; - int nparam_harmonic = fcs->ndup[0].size(); + int nparam_harmonic = fcs->nequiv[0].size(); std::string str_tmp; @@ -391,7 +412,7 @@ void System::load_reference_system() ifs_fc2 >> nat_s >> natmin_ref >> ntran_ref; - if (natmin_ref != symmetry->natmin) { + if (natmin_ref != symmetry->nat_prim) { error->exit("load_reference_system", "The number of atoms in the primitive cell is not consistent"); } @@ -506,10 +527,10 @@ void System::load_reference_system() memory->allocate(ind, 2); list_found.clear(); - for (std::vector::iterator p = fcs->fc_set[0].begin(); - p != fcs->fc_set[0].end(); ++p) { + for (std::vector::iterator p = fcs->fc_table[0].begin(); + p != fcs->fc_table[0].end(); ++p) { for (i = 0; i < 2; ++i) ind[i] = (*p).elems[i]; - list_found.insert(FcProperty(2, (*p).coef, ind, (*p).mother)); + list_found.insert(FcProperty(2, (*p).sign, ind, (*p).mother)); } for (i = 0; i < nparam_harmonic; ++i) { diff --git a/alm/system.h b/alm/system.h index a358760b..450ea699 100644 --- a/alm/system.h +++ b/alm/system.h @@ -68,4 +68,3 @@ namespace ALM_NS void setup_atomic_class(int *); }; } - diff --git a/alm/timer.h b/alm/timer.h index 760923c3..1c963a8d 100644 --- a/alm/timer.h +++ b/alm/timer.h @@ -42,4 +42,3 @@ namespace ALM_NS #endif }; } - diff --git a/alm/writes.cpp b/alm/writes.cpp index d384ec3e..f84e4271 100644 --- a/alm/writes.cpp +++ b/alm/writes.cpp @@ -145,11 +145,11 @@ void Writes::write_force_constants() m = 0; - if (fcs->ndup[order].size() > 0) { + if (fcs->nequiv[order].size() > 0) { ofs_fcs << std::endl << std::setw(6) << str_fcs[order] << std::endl; - for (ui = 0; ui < fcs->ndup[order].size(); ++ui) { + for (ui = 0; ui < fcs->nequiv[order].size(); ++ui) { ofs_fcs << std::setw(8) << k + 1 << std::setw(8) << ui + 1 << std::setw(18) << std::setprecision(7) @@ -157,9 +157,9 @@ void Writes::write_force_constants() atom_tmp.clear(); for (l = 1; l < order + 2; ++l) { - atom_tmp.push_back(fcs->fc_set[order][m].elems[l] / 3); + atom_tmp.push_back(fcs->fc_table[order][m].elems[l] / 3); } - j = symmetry->map_s2p[fcs->fc_set[order][m].elems[0] / 3].atom_num; + j = symmetry->map_s2p[fcs->fc_table[order][m].elems[0] / 3].atom_num; std::sort(atom_tmp.begin(), atom_tmp.end()); iter_cluster = interaction->mindist_cluster[order][j].find( @@ -181,12 +181,12 @@ void Writes::write_force_constants() for (l = 0; l < order + 2; ++l) { ofs_fcs << std::setw(7) - << fcs->easyvizint(fcs->fc_set[order][m].elems[l]); + << fcs->easyvizint(fcs->fc_table[order][m].elems[l]); } ofs_fcs << std::setw(12) << std::setprecision(3) << std::fixed << distmax << std::endl; - m += fcs->ndup[order][ui]; + m += fcs->nequiv[order][ui]; ++k; } } @@ -198,7 +198,7 @@ void Writes::write_force_constants() ofs_fcs << " -------------- Constraints from crystal symmetry --------------" << std::endl << std::endl;; for (order = 0; order < maxorder; ++order) { - int nparam = fcs->ndup[order].size(); + int nparam = fcs->nequiv[order].size(); for (std::vector::iterator p = constraint->const_symmetry[order].begin(); @@ -238,24 +238,24 @@ void Writes::write_force_constants() id = 0; - if (fcs->ndup[order].size() > 0) { + if (fcs->nequiv[order].size() > 0) { ofs_fcs << std::endl << std::setw(6) << str_fcs[order] << std::endl; - for (unsigned int iuniq = 0; iuniq < fcs->ndup[order].size(); ++iuniq) { + for (unsigned int iuniq = 0; iuniq < fcs->nequiv[order].size(); ++iuniq) { str_tmp = " # FC" + boost::lexical_cast(order + 2) + "_"; str_tmp += boost::lexical_cast(iuniq + 1); - ofs_fcs << str_tmp << std::setw(5) << fcs->ndup[order][iuniq] + ofs_fcs << str_tmp << std::setw(5) << fcs->nequiv[order][iuniq] << std::setw(16) << std::scientific << std::setprecision(7) << fitting->params[ip] << std::endl; - for (j = 0; j < fcs->ndup[order][iuniq]; ++j) { + for (j = 0; j < fcs->nequiv[order][iuniq]; ++j) { ofs_fcs << std::setw(5) << j + 1 << std::setw(12) - << std::setprecision(5) << std::fixed << fcs->fc_set[order][id].coef; + << std::setprecision(5) << std::fixed << fcs->fc_table[order][id].sign; for (k = 0; k < order + 2; ++k) { ofs_fcs << std::setw(6) - << fcs->easyvizint(fcs->fc_set[order][id].elems[k]); + << fcs->easyvizint(fcs->fc_table[order][id].elems[k]); } ofs_fcs << std::endl; ++id; @@ -280,8 +280,8 @@ void Writes::write_displacement_pattern() int counter; std::ofstream ofs_pattern; - - std::cout << " Suggested displacement patterns are printed in the following files: " << std::endl; + std::cout << std::endl; + std::cout << " Displacement patterns are saved in the following files: " << std::endl; for (order = 0; order < maxorder; ++order) { ofs_pattern.open(files->file_disp_pattern[order].c_str(), std::ios::out); @@ -293,7 +293,7 @@ void Writes::write_displacement_pattern() ofs_pattern << "Basis : " << displace->disp_basis[0] << std::endl; - for (std::vector::iterator it = displace->pattern_all[order].begin(); + for (auto it = displace->pattern_all[order].begin(); it != displace->pattern_all[order].end(); ++it) { AtomWithDirection entry = *it; @@ -331,7 +331,7 @@ void Writes::write_misc_xml() } system_structure.nat = system->nat; - system_structure.natmin = symmetry->natmin; + system_structure.natmin = symmetry->nat_prim; system_structure.ntran = symmetry->ntran; system_structure.nspecies = system->nkd; @@ -420,7 +420,7 @@ void Writes::write_misc_xml() pt.put("Data.ForceConstants", ""); str_tmp.clear(); - pt.put("Data.ForceConstants.HarmonicUnique.NFC2", fcs->ndup[0].size()); + pt.put("Data.ForceConstants.HarmonicUnique.NFC2", fcs->nequiv[0].size()); int ihead = 0; int k = 0; @@ -429,21 +429,21 @@ void Writes::write_misc_xml() memory->allocate(pair_tmp, nelem); - for (unsigned int ui = 0; ui < fcs->ndup[0].size(); ++ui) { + for (unsigned int ui = 0; ui < fcs->nequiv[0].size(); ++ui) { for (i = 0; i < 2; ++i) { - pair_tmp[i] = fcs->fc_set[0][ihead].elems[i] / 3; + pair_tmp[i] = fcs->fc_table[0][ihead].elems[i] / 3; } j = symmetry->map_s2p[pair_tmp[0]].atom_num; ptree &child = pt.add("Data.ForceConstants.HarmonicUnique.FC2", double2string(fitting->params[k])); child.put(".pairs", - boost::lexical_cast(fcs->fc_set[0][ihead].elems[0]) - + " " + boost::lexical_cast(fcs->fc_set[0][ihead].elems[1])); + boost::lexical_cast(fcs->fc_table[0][ihead].elems[0]) + + " " + boost::lexical_cast(fcs->fc_table[0][ihead].elems[1])); child.put(".multiplicity", interaction->mindist_pairs[pair_tmp[0]][pair_tmp[1]].size()); - ihead += fcs->ndup[0][ui]; + ihead += fcs->nequiv[0][ui]; ++k; } ihead = 0; @@ -455,11 +455,11 @@ void Writes::write_misc_xml() if (interaction->maxorder > 1) { - pt.put("Data.ForceConstants.CubicUnique.NFC3", fcs->ndup[1].size()); + pt.put("Data.ForceConstants.CubicUnique.NFC3", fcs->nequiv[1].size()); - for (unsigned int ui = 0; ui < fcs->ndup[1].size(); ++ui) { + for (unsigned int ui = 0; ui < fcs->nequiv[1].size(); ++ui) { for (i = 0; i < 3; ++i) { - pair_tmp[i] = fcs->fc_set[1][ihead].elems[i] / 3; + pair_tmp[i] = fcs->fc_table[1][ihead].elems[i] / 3; } j = symmetry->map_s2p[pair_tmp[0]].atom_num; @@ -481,21 +481,21 @@ void Writes::write_misc_xml() ptree &child = pt.add("Data.ForceConstants.CubicUnique.FC3", double2string(fitting->params[k])); child.put(".pairs", - boost::lexical_cast(fcs->fc_set[1][ihead].elems[0]) - + " " + boost::lexical_cast(fcs->fc_set[1][ihead].elems[1]) - + " " + boost::lexical_cast(fcs->fc_set[1][ihead].elems[2])); + boost::lexical_cast(fcs->fc_table[1][ihead].elems[0]) + + " " + boost::lexical_cast(fcs->fc_table[1][ihead].elems[1]) + + " " + boost::lexical_cast(fcs->fc_table[1][ihead].elems[2])); child.put(".multiplicity", multiplicity); - ihead += fcs->ndup[1][ui]; + ihead += fcs->nequiv[1][ui]; ++k; } } int ip, ishift; - std::sort(fcs->fc_set[0].begin(), fcs->fc_set[0].end()); + std::sort(fcs->fc_table[0].begin(), fcs->fc_table[0].end()); - for (std::vector::iterator it = fcs->fc_set[0].begin(); - it != fcs->fc_set[0].end(); ++it) { + for (std::vector::iterator it = fcs->fc_table[0].begin(); + it != fcs->fc_table[0].end(); ++it) { FcProperty fctmp = *it; ip = fctmp.mother; @@ -506,7 +506,7 @@ void Writes::write_misc_xml() for (std::vector::iterator it2 = interaction->mindist_pairs[pair_tmp[0]][pair_tmp[1]].begin(); it2 != interaction->mindist_pairs[pair_tmp[0]][pair_tmp[1]].end(); ++it2) { ptree &child = pt.add("Data.ForceConstants.HARMONIC.FC2", - double2string(fitting->params[ip] * fctmp.coef + double2string(fitting->params[ip] * fctmp.sign / static_cast(interaction->mindist_pairs[pair_tmp[0]][pair_tmp[1]].size()))); child.put(".pair1", boost::lexical_cast(j + 1) @@ -517,7 +517,7 @@ void Writes::write_misc_xml() } } - ishift = fcs->ndup[0].size(); + ishift = fcs->nequiv[0].size(); // Print anharmonic force constants to the xml file. @@ -527,10 +527,10 @@ void Writes::write_misc_xml() std::string elementname; for (order = 1; order < interaction->maxorder; ++order) { - std::sort(fcs->fc_set[order].begin(), fcs->fc_set[order].end()); + std::sort(fcs->fc_table[order].begin(), fcs->fc_table[order].end()); - for (std::vector::iterator it = fcs->fc_set[order].begin(); - it != fcs->fc_set[order].end(); ++it) { + for (std::vector::iterator it = fcs->fc_table[order].begin(); + it != fcs->fc_table[order].end(); ++it) { FcProperty fctmp = *it; ip = fctmp.mother + ishift; @@ -561,7 +561,7 @@ void Writes::write_misc_xml() std::vector cell_now = (*iter_cluster).cell[imult]; ptree &child = pt.add(elementname, - double2string(fitting->params[ip] * fctmp.coef + double2string(fitting->params[ip] * fctmp.sign / static_cast(multiplicity))); child.put(".pair1", boost::lexical_cast(j + 1) @@ -578,7 +578,7 @@ void Writes::write_misc_xml() error->exit("write_misc_xml", "This cannot happen."); } } - ishift += fcs->ndup[order].size(); + ishift += fcs->nequiv[order].size(); } using namespace boost::property_tree::xml_parser; @@ -617,8 +617,8 @@ void Writes::write_hessian() } } - for (std::vector::iterator it = fcs->fc_set[0].begin(); - it != fcs->fc_set[0].end(); ++it) { + for (std::vector::iterator it = fcs->fc_table[0].begin(); + it != fcs->fc_table[0].end(); ++it) { FcProperty fctmp = *it; ip = fctmp.mother; @@ -628,7 +628,7 @@ void Writes::write_hessian() pair_tran[i] = symmetry->map_sym[pair_tmp[i]][symmetry->symnum_tran[itran]]; } hessian[3 * pair_tran[0] + fctmp.elems[0] % 3][3 * pair_tran[1] + fctmp.elems[1] % 3] - = fitting->params[ip] * fctmp.coef; + = fitting->params[ip] * fctmp.sign; } } @@ -658,8 +658,8 @@ void Writes::write_hessian() ofs_fc2.open(file_fc2.c_str(), std::ios::out); ofs_fc2 << " # iat, icrd, jat, icrd, icell, relvec, fc2" << std::endl; double vec[3]; - for (std::vector::iterator it = fcs->fc_set[0].begin(); - it != fcs->fc_set[0].end(); ++it) { + for (std::vector::iterator it = fcs->fc_table[0].begin(); + it != fcs->fc_table[0].end(); ++it) { FcProperty fctmp = *it; ip = fctmp.mother; @@ -684,7 +684,7 @@ void Writes::write_hessian() ofs_fc2 << std::setw(15) << vec[2]; ofs_fc2 << std::setw(15) - << fitting->params[ip] * fctmp.coef / static_cast(multiplicity); + << fitting->params[ip] * fctmp.sign / static_cast(multiplicity); ofs_fc2 << std::endl; } } diff --git a/alm/writes.h b/alm/writes.h index 0e48f7b8..7225e6bc 100644 --- a/alm/writes.h +++ b/alm/writes.h @@ -24,11 +24,15 @@ namespace ALM_NS int kind; int atom, tran; - AtomProperty() {}; + AtomProperty() + { + }; AtomProperty(const AtomProperty &other) : x(other.x), y(other.y), z(other.z), - kind(other.kind), atom(other.atom), tran(other.tran) {}; + kind(other.kind), atom(other.atom), tran(other.tran) + { + }; AtomProperty(const double *pos, const int kind_in, @@ -52,7 +56,9 @@ namespace ALM_NS int nat, natmin, ntran; int nspecies; - SystemInfo() {}; + SystemInfo() + { + }; }; class Writes: protected Pointers diff --git a/anphon/anphon.vcxproj b/anphon/anphon.vcxproj index 4390d6a1..9661b3f3 100644 --- a/anphon/anphon.vcxproj +++ b/anphon/anphon.vcxproj @@ -43,7 +43,7 @@ $(UniversalCRT_LibraryPath_x86);$(MSMPI_LIB32);C:\fftw3\lib;$(VCInstallDir)lib;$(VCInstallDir)atlmfc\lib;$(WindowsSDK_LibraryPath_x86) - $(UniversalCRT_IncludePath);$(MSMPI_INC);$(MSMPI_INC)x86;C:\fftw3\include;C:\boost\boost_1_60_0;$(VCInstallDir)include;$(VCInstallDir)atlmfc\include;$(WindowsSDK_IncludePath) + C:\eigen_c++;$(UniversalCRT_IncludePath);$(MSMPI_INC);$(MSMPI_INC)x86;C:\fftw3\include;C:\boost\boost_1_60_0;$(VCInstallDir)include;$(VCInstallDir)atlmfc\include;$(WindowsSDK_IncludePath) $(UniversalCRT_LibraryPath_x86);$(MSMPI_LIB32);C:\fftw3\lib;$(VCInstallDir)lib;$(VCInstallDir)atlmfc\lib;$(WindowsSDK_LibraryPath_x86); diff --git a/anphon/dynamical.cpp b/anphon/dynamical.cpp index 94c35e2a..4699634f 100644 --- a/anphon/dynamical.cpp +++ b/anphon/dynamical.cpp @@ -28,11 +28,15 @@ #include "phonon_dos.h" #include "gruneisen.h" #include "ewald.h" +#include + using namespace PHON_NS; Dynamical::Dynamical(PHON *phon): Pointers(phon) { + index_bconnect = nullptr; + symmetrize_borncharge = 0; } Dynamical::~Dynamical() @@ -113,11 +117,12 @@ 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); - if (mympi->my_rank == 0) load_born(); + if (mympi->my_rank == 0) load_born(symmetrize_borncharge); MPI_Bcast(&dielec[0][0], 9, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(&borncharge[0][0][0], 9 * system->natmin, MPI_DOUBLE, 0, MPI_COMM_WORLD); @@ -147,6 +152,10 @@ void PHON_NS::Dynamical::finish_dynamical() if (nonanalytic) { memory->deallocate(borncharge); } + + if (index_bconnect) { + memory->deallocate(index_bconnect); + } } @@ -740,6 +749,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; } @@ -804,7 +818,7 @@ void Dynamical::modify_eigenvectors() } -void Dynamical::load_born() +void Dynamical::load_born(const unsigned int flag_symmborn) { // Read the dielectric tensor and born effective charges from file_born @@ -879,7 +893,7 @@ void Dynamical::load_born() if (res > eps10) { std::cout << std::endl; std::cout << " WARNING: Born effective charges do not satisfy the acoustic sum rule." << std::endl; - std::cout << " The born effective charges are modified to follow the ASR." << std::endl; + std::cout << " The born effective charges are modified to satisfy the ASR." << std::endl; for (i = 0; i < system->natmin; ++i) { for (j = 0; j < 3; ++j) { @@ -890,91 +904,94 @@ void Dynamical::load_born() } } - // Symmetrize Born effective charges. Necessary to avoid the violation of ASR - // particularly for NONANALYTIC=3 (Ewald summation). + if (flag_symmborn) { - int isym, iat, iat_sym; - int m; - double ***born_sym; - double rot[3][3]; + // Symmetrize Born effective charges. Necessary to avoid the violation of ASR + // particularly for NONANALYTIC=3 (Ewald summation). - memory->allocate(born_sym, system->natmin, 3, 3); + int isym, iat, iat_sym; + int m; + double ***born_sym; + double rot[3][3]; - for (iat = 0; iat < system->natmin; ++iat) { - for (i = 0; i < 3; ++i) { - for (j = 0; j < 3; ++j) { - born_sym[iat][i][j] = 0.0; - } - } - } + memory->allocate(born_sym, system->natmin, 3, 3); - for (isym = 0; isym < symmetry->SymmListWithMap.size(); ++isym) { - for (i = 0; i < 3; ++i) { - for (j = 0; j < 3; ++j) { - rot[i][j] = symmetry->SymmListWithMap[isym].rot[3 * i + j]; + for (iat = 0; iat < system->natmin; ++iat) { + for (i = 0; i < 3; ++i) { + for (j = 0; j < 3; ++j) { + born_sym[iat][i][j] = 0.0; + } } } - for (iat = 0; iat < system->natmin; ++iat) { - iat_sym = symmetry->SymmListWithMap[isym].mapping[iat]; - + for (isym = 0; isym < symmetry->SymmListWithMap.size(); ++isym) { for (i = 0; i < 3; ++i) { for (j = 0; j < 3; ++j) { - for (k = 0; k < 3; ++k) { - for (m = 0; m < 3; ++m) { - born_sym[iat][i][j] += rot[i][k] * rot[j][m] * borncharge[iat_sym][k][m]; + rot[i][j] = symmetry->SymmListWithMap[isym].rot[3 * i + j]; + } + } + + for (iat = 0; iat < system->natmin; ++iat) { + iat_sym = symmetry->SymmListWithMap[isym].mapping[iat]; + + for (i = 0; i < 3; ++i) { + for (j = 0; j < 3; ++j) { + for (k = 0; k < 3; ++k) { + for (m = 0; m < 3; ++m) { + born_sym[iat_sym][i][j] += rot[i][k] * rot[j][m] * borncharge[iat][k][m]; + } } } } } } - } - for (iat = 0; iat < system->natmin; ++iat) { - for (i = 0; i < 3; ++i) { - for (j = 0; j < 3; ++j) { - born_sym[iat][i][j] /= static_cast(symmetry->SymmListWithMap.size()); + for (iat = 0; iat < system->natmin; ++iat) { + for (i = 0; i < 3; ++i) { + for (j = 0; j < 3; ++j) { + born_sym[iat][i][j] /= static_cast(symmetry->SymmListWithMap.size()); + } } } - } - // Check if the Born effective charges given by the users satisfy the symmetry. + // Check if the Born effective charges given by the users satisfy the symmetry. - double diff_sym = 0.0; - for (iat = 0; iat < system->natmin; ++iat) { - for (i = 0; i < 3; ++i) { - for (j = 0; j < 3; ++j) { - diff_sym = std::max(res, std::abs(borncharge[iat][i][j] - born_sym[iat][i][j])); + double diff_sym = 0.0; + for (iat = 0; iat < system->natmin; ++iat) { + for (i = 0; i < 3; ++i) { + for (j = 0; j < 3; ++j) { + diff_sym = std::max(diff_sym, std::abs(borncharge[iat][i][j] - born_sym[iat][i][j])); + } } } - } - if (diff_sym > 0.5) { - std::cout << std::endl; - std::cout << " WARNING: Born effective charges are inconsistent with the crystal symmetry." << std::endl; - } + if (diff_sym > 0.5) { + std::cout << std::endl; + std::cout << " WARNING: Born effective charges are inconsistent with the crystal symmetry." << std::endl; + } - for (iat = 0; iat < system->natmin; ++iat) { - for (i = 0; i < 3; ++i) { - for (j = 0; j < 3; ++j) { - borncharge[iat][i][j] = born_sym[iat][i][j]; + for (iat = 0; iat < system->natmin; ++iat) { + for (i = 0; i < 3; ++i) { + for (j = 0; j < 3; ++j) { + borncharge[iat][i][j] = born_sym[iat][i][j]; + } } } - } - memory->deallocate(born_sym); + 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; + 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]; + for (j = 0; j < 3; ++j) { + for (k = 0; k < 3; ++k) { + std::cout << std::setw(15) << borncharge[i][j][k]; + } + std::cout << std::endl; } - std::cout << std::endl; } } } @@ -1054,3 +1071,90 @@ void Dynamical::calc_atomic_participation_ratio(std::complex *evec, doub for (iat = 0; iat < natmin; ++iat) ret[iat] /= std::sqrt(static_cast(natmin) * sum); } + + +void Dynamical::connect_band_by_eigen_similarity(std::complex ***evec, + int **index_sorted) +{ + int ik, is, js, ks; + unsigned int nk = kpoint->nk; + unsigned int ns = neval; + int loc; + std::vector index; + std::complex **evec_tmp; + std::vector> abs_similarity; + std::complex dprod; + std::vector 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(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) { + + // 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); +} diff --git a/anphon/dynamical.h b/anphon/dynamical.h index 9d17fc6f..9d84b345 100644 --- a/anphon/dynamical.h +++ b/anphon/dynamical.h @@ -51,13 +51,16 @@ namespace PHON_NS unsigned int neval; bool eigenvectors; bool print_eigenvectors; + unsigned int symmetrize_borncharge; 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 ***evec_phonon; double dielec[3][3]; double ***borncharge; @@ -96,11 +99,12 @@ namespace PHON_NS private: - void load_born(); + void load_born(const unsigned int); void prepare_mindist_list(std::vector **); void calc_atomic_participation_ratio(std::complex *, double *); double distance(double *, double *); + void connect_band_by_eigen_similarity(std::complex ***, int **); // void calc_analytic_k(double *, double ****, std::complex **); // void modify_eigenvectors_sym(); diff --git a/anphon/parsephon.cpp b/anphon/parsephon.cpp index e8caccd9..59bbcfa3 100644 --- a/anphon/parsephon.cpp +++ b/anphon/parsephon.cpp @@ -120,6 +120,8 @@ void Input::parse_general_vars() bool selenergy_offdiagonal; bool update_fc2; bool classical; + unsigned int bornsym; + unsigned int band_connection; struct stat st; std::string prefix, mode, fcsinfo, fc2info; @@ -128,7 +130,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 BORNSYM"; std::string str_no_defaults = "PREFIX MODE FCSXML NKD KD MASS"; std::vector no_defaults; std::vector kdname_v, masskd_v; @@ -207,6 +209,8 @@ void Input::parse_general_vars() sym_time_reversal = false; use_triplet_symmetry = true; classical = false; + band_connection = 0; + bornsym = 0; prec_ewald = 1.0e-12; @@ -253,8 +257,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); + assign_val(bornsym, "BORNSYM", general_var_dict); + + if (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); @@ -312,8 +321,10 @@ void Input::parse_general_vars() dynamical->nonanalytic = nonanalytic; dynamical->na_sigma = na_sigma; + dynamical->symmetrize_borncharge = bornsym; writes->nbands = nbands; dynamical->file_born = borninfo; + dynamical->band_connection = band_connection; integration->epsilon = epsilon; fcs_phonon->file_fcs = fcsinfo; if (!fc2info.empty()) { diff --git a/anphon/pointers.h b/anphon/pointers.h index ac2946ca..6b681003 100644 --- a/anphon/pointers.h +++ b/anphon/pointers.h @@ -40,9 +40,13 @@ namespace PHON_NS isotope(ptr->isotope), scph(ptr->scph), ewald(ptr->ewald), - timer(ptr->timer) {} + timer(ptr->timer) + { + } - virtual ~Pointers() {} + virtual ~Pointers() + { + } protected: PHON *phon; @@ -70,4 +74,3 @@ namespace PHON_NS Timer *&timer; }; } - diff --git a/anphon/relaxation.cpp b/anphon/relaxation.cpp index b1ae48e4..43be072d 100644 --- a/anphon/relaxation.cpp +++ b/anphon/relaxation.cpp @@ -1208,7 +1208,7 @@ void Relaxation::calc_frequency_resolved_final_state(const unsigned int N, 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) @@ -2542,7 +2542,7 @@ void Relaxation::print_momentum_resolved_final_state(const unsigned int NT, n1 = f1 + f2 + 1.0; n2 = f1 - f2; } - + if (selection_type == 0) { gamma_k[k][iT] += V3norm * n1; @@ -3120,7 +3120,7 @@ void Relaxation::calc_self3omega_tetrahedron(const double 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]); diff --git a/anphon/scph.cpp b/anphon/scph.cpp index 0b9bb89d..1ee007fa 100644 --- a/anphon/scph.cpp +++ b/anphon/scph.cpp @@ -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); } diff --git a/anphon/selfenergy.cpp b/anphon/selfenergy.cpp index c1d161e5..b1cf66a1 100644 --- a/anphon/selfenergy.cpp +++ b/anphon/selfenergy.cpp @@ -127,7 +127,7 @@ void Selfenergy::selfenergy_tadpole(const unsigned int N, } else { n2 = thermodynamics->fB(omega2, T_tmp); ret_mpi[i] += v3_tmp2 * (2.0 * n2 + 1.0); - } + } } } } diff --git a/anphon/thermodynamics.cpp b/anphon/thermodynamics.cpp index f9f94d54..cb57a4cb 100644 --- a/anphon/thermodynamics.cpp +++ b/anphon/thermodynamics.cpp @@ -273,7 +273,7 @@ double Thermodynamics::free_energy(const double T) x = omega / (T * T_to_Ryd); ret += std::log(x); } - + return T * T_to_Ryd * ret / static_cast(nk); } else { @@ -297,7 +297,6 @@ double Thermodynamics::free_energy(const double T) return T * T_to_Ryd * ret / static_cast(nk); } - } double Thermodynamics::disp2_avg(const double T, @@ -325,7 +324,7 @@ double Thermodynamics::disp2_avg(const double T, if (omega < eps8) continue; ret += real(dynamical->evec_phonon[ik][is][ns1] - * std::conj(dynamical->evec_phonon[ik][is][ns2])) + * std::conj(dynamical->evec_phonon[ik][is][ns2])) * T * T_to_Ryd / (omega * omega); } } else { @@ -340,7 +339,7 @@ double Thermodynamics::disp2_avg(const double T, if (omega < eps8) continue; ret += real(dynamical->evec_phonon[ik][is][ns1] - * std::conj(dynamical->evec_phonon[ik][is][ns2])) + * std::conj(dynamical->evec_phonon[ik][is][ns2])) * (fB(omega, T) + 0.5) / omega; } } diff --git a/anphon/write_phonons.cpp b/anphon/write_phonons.cpp index b21d633b..17a192bd 100644 --- a/anphon/write_phonons.cpp +++ b/anphon/write_phonons.cpp @@ -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") { @@ -281,7 +282,7 @@ void Writes::setup_result_io() } if (static_cast(is_classical) != thermodynamics->classical) { error->warn("setup_result_io", - "CLASSICAL val is not consistent"); + "CLASSICAL val is not consistent"); } found_tag = false; @@ -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() diff --git a/docs/source/conf.py b/docs/source/conf.py index e4076b56..67eb28ed 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -62,7 +62,7 @@ # The short X.Y version. version = '1.0' # The full version, including alpha/beta/rc tags. -release = '1.0.0' +release = '1.0.2' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/docs/source/formalism/formalism_alm.rst b/docs/source/formalism/formalism_alm.rst index d936d6d9..aec8c3a9 100644 --- a/docs/source/formalism/formalism_alm.rst +++ b/docs/source/formalism/formalism_alm.rst @@ -45,14 +45,14 @@ The are several relationships between IFCs which may be used to reduce the numbe .. math:: :label: ifcsym1 - \sum_{\nu_{1},\dots,\nu_{n}}\Phi_{\nu_{1}\dots\nu_{n}}(L_{1}K_{1};\dots;L_{n}K_{n}) O_{\mu_{1}\nu_{1}}\cdots O_{\mu_{n}\nu_{n}} = \Phi_{\mu_{1}\dots\mu_{n}}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n}), + \sum_{\nu_{1},\dots,\nu_{n}}\Phi_{\nu_{1}\dots\nu_{n}}(L_{1}K_{1};\dots;L_{n}K_{n}) O_{\nu_{1}\mu_{1}}\cdots O_{\nu_{n}\mu_{n}} = \Phi_{\mu_{1}\dots\mu_{n}}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n}), where :math:`O` is the rotational matrix of the symmetry operation. Let :math:`N_{s}` be the number of symmetry operations, there are :math:`N_{s}` relationships between IFCs which may be used to find independent IFCs. .. Note:: - In the current implementation of *alm*, independent IFCs are searched in Cartesian coordinate where the matrix element :math:`O_{\mu\nu}` is 0 or :math:`\pm1` in all symmetry operations except for those of **hexagonal** (trigonal) lattice. Also, except for hexagonal (trigonal) systems, the product :math:`O_{\mu_{1}\nu_{1}}\cdots O_{\mu_{n}\nu_{n}}` in the left hand side of equation :eq:`ifcsym1` becomes non-zero only for a specific pair of :math:`\{\nu\}` (and becomes 0 for all other :math:`\{\nu\}`\ s). Therefore, let :math:`\{\nu^{\prime}\}` be such a pair of :math:`\{\nu\}`, the equation :eq:`ifcsym1` can be reduced to + In the current implementation of *alm*, independent IFCs are searched in Cartesian coordinate where the matrix element :math:`O_{\mu\nu}` is 0 or :math:`\pm1` in all symmetry operations except for those of **hexagonal** (trigonal) lattice. Also, except for hexagonal (trigonal) systems, the product :math:`O_{\nu_{1}\mu_{1}}\cdots O_{\nu_{n}\mu_{n}}` in the left hand side of equation :eq:`ifcsym1` becomes non-zero only for a specific pair of :math:`\{\nu\}` (and becomes 0 for all other :math:`\{\nu\}`\ s). Therefore, let :math:`\{\nu^{\prime}\}` be such a pair of :math:`\{\nu\}`, the equation :eq:`ifcsym1` can be reduced to .. math:: :label: ifcsym2 diff --git a/docs/source/input/inputanphon.rst b/docs/source/input/inputanphon.rst index f5488398..ebd938cf 100644 --- a/docs/source/input/inputanphon.rst +++ b/docs/source/input/inputanphon.rst @@ -114,7 +114,6 @@ List of input variables === ======================================================= 0 Symmetry operations won’t be saved in “SYMM_INFO_PRIM” - 1 Symmetry operations will be saved in “SYMM_INFO_PRIM” === ======================================================= @@ -157,6 +156,18 @@ List of input variables ```` +* BORNSYM-tag = 0 | 1 + + === ================================================================= + 0 Do not symmetrize Born effective charges + 1 Symmetrize Born effective charges by using point group symmetry + === ================================================================= + + :Default: 0 + :Type: Integer + +```` + * TMIN, TMAX, DT-tags : Temperature range and its stride in units of Kelvin :Default: ``TMIN = 0``, ``TMAX = 1000``, ``DT = 10`` @@ -193,6 +204,25 @@ List of input variables ```` +* BCONNECT-tag = 0 | 1 | 2 + + === =================================================================================== + 0 | Phonon band is saved without change (sorted in order of energy) + + 1 | Phonon band is connected by using the similarity of eigenvectors. + + 2 | Same as ``BCONNECT=1``. In addition, information of the connectivity is + | saved as ``PREFIX.connection``. + === =================================================================================== + + :Default: 0 + :Type: Integer + :Description: The algorithm for connecting a band structure is described here_. + + .. _here : https://www.slideshare.net/TakeshiNishimatsu/two-efficient-algorithms-for-drawing-accurate-and-beautiful-phonon-dispersion + +```` + * CLASSICAL-tag = 0 | 1 === ======================================================= diff --git a/include/mathfunctions.h b/include/mathfunctions.h index 786f410f..26c0ad6e 100644 --- a/include/mathfunctions.h +++ b/include/mathfunctions.h @@ -14,7 +14,7 @@ #include template -inline void matmul3(T ret[3][3], T amat[3][3], T bmat[3][3]) { +inline void matmul3(T ret[3][3], const T amat[3][3], const T bmat[3][3]) { int i, j, k; T ret_tmp[3][3]; diff --git a/include/version.h b/include/version.h index cae9d288..866d08e6 100644 --- a/include/version.h +++ b/include/version.h @@ -12,4 +12,4 @@ #include -static const std::string ALAMODE_VERSION = "1.0.0"; +static const std::string ALAMODE_VERSION = "1.0.2"; diff --git a/tools/analyze_phonons.cpp b/tools/analyze_phonons.cpp index 917878d0..e9352dc9 100644 --- a/tools/analyze_phonons.cpp +++ b/tools/analyze_phonons.cpp @@ -563,7 +563,7 @@ void calc_kappa_cumulative(double max_length, double delta_length, int itemp) } else { c_tmp = Cv(omega[ik][is], temp[itemp]); } - + vel_tmp = pow(vel[ik][is][0][0], 2) + pow(vel[ik][is][0][1], 2) + pow(vel[ik][is][0][2], 2); @@ -782,7 +782,7 @@ void calc_kappa_boundary2(double max_length, double delta_length, int itemp, int c_tmp = k_Boltzmann; } else { c_tmp = Cv(omega[ik][is], temp[itemp]); - } + } for (i = 0; i < nsame; ++i) { diff --git a/tools/extract.py b/tools/extract.py index a2009262..1ea29bb3 100755 --- a/tools/extract.py +++ b/tools/extract.py @@ -124,7 +124,7 @@ def print_displacements_VASP(xml_files, for search_target in xml_files: x = get_coordinate_VASP(search_target, nat) - ndata = len(x) / (3 * nat) + ndata = len(x) // (3 * nat) x = np.reshape(x, (ndata, nat, 3)) for idata in range(ndata): @@ -179,7 +179,7 @@ def print_atomicforces_VASP(xml_files, for search_target in xml_files: data = get_atomicforces_VASP(search_target) - ndata = len(data) / (3 * nat) + ndata = len(data) // (3 * nat) data = np.reshape(data, (ndata, nat, 3)) for idata in range(ndata):