Skip to content

Commit

Permalink
For fc_cart and get_fc_cart, the zero elements should be kept and
Browse files Browse the repository at this point in the history 
returned because it will be used to fix force constants. In some
cases, the force constants should be fixed to zero due to symmetry
constraints. To properly set the constraint matrix in such cases,
the get_fc method in the python API should also return zero elements.
  • Loading branch information
@ttadano
ttadano committed Mar 30, 2023
1 parent bdc19d8 commit 96b32e5
Show file tree
Hide file tree
Showing 4 changed files with 38 additions and 53 deletions.
47 changes: 14 additions & 33 deletions src/constraint.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2743,30 +2743,21 @@ void Constraint::generate_fix_constraint(const Symmetry *symmetry,
fcs->translate_forceconstant_index_to_centercell(symmetry,
intpair_to_fix);

// std::cout << "intpair_fix.size() = " << intpair_to_fix.size() << '\n';
// for (auto ifix = 0; ifix < intpair_to_fix.size(); ++ifix) {
// for (auto it: intpair_to_fix[ifix]) {
// std::cout << std::setw(5) << it;
// }
// std::cout << std::setw(15) << values_fix_fc2[ifix] << '\n';
// }

std::vector<ForceConstantTable> fc_fix_table;
std::set<ForceConstantTable> fc_fix_table;

const auto nfcs = intpair_to_fix.size();

for (auto i = 0; i < nfcs; ++i) {
fc_fix_table.emplace_back(values_fix_fc2[i],
intpair_to_fix[i]);
fc_fix_table.insert(ForceConstantTable(values_fix_fc2[i],
intpair_to_fix[i]));
}
std::sort(fc_fix_table.begin(), fc_fix_table.end());

size_t ihead = 0;

std::vector<int> index_tmp;
double sign;
size_t mother;
std::vector<ForceConstantTable>::iterator it_found;
std::set<ForceConstantTable>::iterator it_found;
bool found_element;

const_fix[order].clear();
Expand All @@ -2783,9 +2774,7 @@ void Constraint::generate_fix_constraint(const Symmetry *symmetry,
for (auto j = 0; j < fcs->get_nequiv()[order][ui]; ++j) {
index_tmp = fcs->get_fc_table()[order][ihead + j].elems;

it_found = std::lower_bound(fc_fix_table.begin(),
fc_fix_table.end(),
ForceConstantTable(0.0, index_tmp));
it_found = fc_fix_table.find(ForceConstantTable(0.0, index_tmp));

if (it_found != fc_fix_table.end()) {
found_element = true;
Expand All @@ -2796,15 +2785,10 @@ void Constraint::generate_fix_constraint(const Symmetry *symmetry,
break;
}
}
//
// std::cout << "found: mother = " << std::setw(5) << mother;
// std::cout << " val = " << std::setw(15) << sign * (*it_found).fc_value;
// std::cout << " sign_mother = " << sign_mother;
// std::cout << '\n';

if (found_element) {
const_fix[order].emplace_back(ConstraintTypeFix(mother,
sign * (*it_found).fc_value));
const_fix[order].emplace_back(mother,
sign * (*it_found).fc_value);
}
ihead += fcs->get_nequiv()[order][ui];
}
Expand All @@ -2820,20 +2804,19 @@ void Constraint::generate_fix_constraint(const Symmetry *symmetry,

const auto nfcs = intpair_to_fix.size();

std::vector<ForceConstantTable> fc_fix_table;
std::set<ForceConstantTable> fc_fix_table;

for (auto i = 0; i < nfcs; ++i) {
fc_fix_table.emplace_back(values_fix_fc3[i],
intpair_to_fix[i]);
fc_fix_table.insert(ForceConstantTable(values_fix_fc3[i],
intpair_to_fix[i]));
}
std::sort(fc_fix_table.begin(), fc_fix_table.end());

size_t ihead = 0;

std::vector<int> index_tmp;
double sign;
size_t mother;
std::vector<ForceConstantTable>::iterator it_found;
std::set<ForceConstantTable>::iterator it_found;
bool found_element;

const_fix[order].clear();
Expand All @@ -2849,9 +2832,7 @@ void Constraint::generate_fix_constraint(const Symmetry *symmetry,

for (auto j = 0; j < fcs->get_nequiv()[order][ui]; ++j) {
index_tmp = fcs->get_fc_table()[order][ihead + j].elems;
it_found = std::lower_bound(fc_fix_table.begin(),
fc_fix_table.end(),
ForceConstantTable(0.0, index_tmp));
it_found = fc_fix_table.find(ForceConstantTable(0.0, index_tmp));

if (it_found != fc_fix_table.end()) {
found_element = true;
Expand All @@ -2864,8 +2845,8 @@ void Constraint::generate_fix_constraint(const Symmetry *symmetry,
}

if (found_element) {
const_fix[order].emplace_back(ConstraintTypeFix(mother,
sign * (*it_found).fc_value));
const_fix[order].emplace_back(mother,
sign * (*it_found).fc_value);
}
ihead += fcs->get_nequiv()[order][ui];
}
Expand Down
36 changes: 16 additions & 20 deletions src/fcs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -903,12 +903,10 @@ void Fcs::set_forceconstant_cartesian(const int maxorder,
fcs_cart += prod_matrix * fc_list_grp[igrp][j];
}

if (std::abs(fcs_cart) >= fc_zero_threshold) {
fc_cart_omp.emplace_back(nelems,
fcs_cart,
&atoms_grp[igrp][0],
xyzcomponent[ixyz]);
}
fc_cart_omp.emplace_back(nelems,
fcs_cart,
&atoms_grp[igrp][0],
xyzcomponent[ixyz]);
}
}

Expand Down Expand Up @@ -939,20 +937,18 @@ void Fcs::set_forceconstant_cartesian(const int maxorder,
coords_now.resize(nelems);

for (const auto &it: fc_table[i]) {
if (std::abs(param_in[it.mother + ishift]) >= fc_zero_threshold) {
elems_permutation = it.elems;
do {
for (j = 0; j < nelems; ++j) {
atoms_now[j] = elems_permutation[j] / 3;
coords_now[j] = elems_permutation[j] % 3;
}
fc_cart[i].emplace_back(nelems,
param_in[it.mother + ishift] * it.sign,
&atoms_now[0],
&coords_now[0]);
} while (std::next_permutation(elems_permutation.begin() + 1,
elems_permutation.end()));
}
elems_permutation = it.elems;
do {
for (j = 0; j < nelems; ++j) {
atoms_now[j] = elems_permutation[j] / 3;
coords_now[j] = elems_permutation[j] % 3;
}
fc_cart[i].emplace_back(nelems,
param_in[it.mother + ishift] * it.sign,
&atoms_now[0],
&coords_now[0]);
} while (std::next_permutation(elems_permutation.begin() + 1,
elems_permutation.end()));
}

ishift += nequiv[i].size();
Expand Down
5 changes: 5 additions & 0 deletions src/fcs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,6 +148,11 @@ class ForceConstantTable {
return std::lexicographical_compare(flattenarray.begin(), flattenarray.end(),
a.flattenarray.begin(), a.flattenarray.end());
}

bool operator==(const ForceConstantTable &a) const
{
return flattenarray == a.flattenarray;
}
};

class Fcs {
Expand Down
3 changes: 3 additions & 0 deletions src/writer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -662,6 +662,8 @@ void Writer::save_fcs_alamode_oldformat(const System *system,

for (const auto &it: fc_cart_harmonic) {

if (std::abs(it.fc_value) < fcs->get_fc_zero_threshold()) continue;

for (k = 0; k < 2; ++k) {
pair_tmp[k] = it.atoms[k];
}
Expand Down Expand Up @@ -717,6 +719,7 @@ void Writer::save_fcs_alamode_oldformat(const System *system,
// and the last (order + 1) elements are sorted in ascending order.

if (!it.is_ascending_order) continue;
if (std::abs(it.fc_value) < fcs->get_fc_zero_threshold()) continue;

for (k = 0; k < order + 2; ++k) {
pair_tmp[k] = it.atoms[k];
Expand Down

0 comments on commit 96b32e5

Please sign in to comment.