Skip to content

perf the compute of check_atom_stru and add openmp #5962

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 17 commits into from
Mar 7, 2025
+156 −140
Merged

perf the compute of check_atom_stru and add openmp #5962

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
0d17a87
first version
A-006 Mar 3, 2025
cfd3778
change add
A-006 Mar 3, 2025
442e95c
add the change
A-006 Mar 3, 2025
215ac56
change
A-006 Mar 3, 2025
471fa50
change gtest
A-006 Mar 4, 2025
3b2e1f2
add change
A-006 Mar 5, 2025
a407f06
change the file
A-006 Mar 5, 2025
d268a3e
add format
A-006 Mar 5, 2025
2ad9f1e
modify back test
A-006 Mar 5, 2025
8e7fa1b
Merge branch 'develop' into ucell27
A-006 Mar 5, 2025
f5adb78
add the change in the check_atom_stru
A-006 Mar 6, 2025
0dde41d
move openmp
A-006 Mar 6, 2025
1101104
add globlv for the question
A-006 Mar 6, 2025
0849de3
Merge branch 'develop' into ucell27
A-006 Mar 6, 2025
e03877c
Merge branch 'develop' into ucell27
A-006 Mar 7, 2025
4e9acb6
delete vector
A-006 Mar 7, 2025
9a91913
update bug
A-006 Mar 7, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
270 changes: 139 additions & 131 deletions source/module_cell/check_atomic_stru.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,157 +1,165 @@
#include "check_atomic_stru.h"

#include "module_base/element_covalent_radius.h"
#include "module_base/timer.h"

void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) {
void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor)
{
// First we calculate all bond length in the structure,
// and compare with the covalent_bond_length,
// if there has bond length is shorter than covalent_bond_length * factor,
// we think this structure is unreasonable.
const double warning_coef = 0.6;
assert(ucell.ntype > 0);
std::stringstream errorlog;
bool all_pass = true;
bool no_warning = true;
for (int it1 = 0; it1 < ucell.ntype; it1++) {
std::string symbol1 = "";
for (char ch: ucell.atoms[it1].label) {
if (std::isalpha(ch)) {
symbol1.push_back(ch);
}
}
// std::string symbol1 = ucell.atoms[it1].label;
double symbol1_covalent_radius;
if (ModuleBase::CovalentRadius.find(symbol1)
!= ModuleBase::CovalentRadius.end()) {
symbol1_covalent_radius = ModuleBase::CovalentRadius.at(symbol1);
} else {
std::stringstream mess;
mess << "Notice: symbol '" << symbol1
<< "' is not an element symbol!!!! ";
mess << "set the covalent radius to be 0." << std::endl;
GlobalV::ofs_running << mess.str();
std::cout << mess.str();
symbol1_covalent_radius = 0.0;
}
std::stringstream errorlog;
errorlog.setf(std::ios_base::fixed, std::ios_base::floatfield);

for (int ia1 = 0; ia1 < ucell.atoms[it1].na; ia1++) {
double x1 = ucell.atoms[it1].taud[ia1].x;
double y1 = ucell.atoms[it1].taud[ia1].y;
double z1 = ucell.atoms[it1].taud[ia1].z;
if (GlobalV::MY_RANK == 0)
{
ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru");

const int ntype = ucell.ntype;
const double lat0 = ucell.lat0;
const double warning_coef = 0.6;
const double max_factor_coef = std::max(warning_coef, factor);

for (int it2 = 0; it2 < ucell.ntype; it2++) {
std::string symbol2 = ucell.atoms[it2].label;
double symbol2_covalent_radius;
if (ModuleBase::CovalentRadius.find(symbol2)
!= ModuleBase::CovalentRadius.end()) {
symbol2_covalent_radius
= ModuleBase::CovalentRadius.at(symbol2);
} else {
symbol2_covalent_radius = 0.0;
std::vector<double> symbol_covalent_radiuss(ntype);
for (int it = 0; it < ntype; it++)
{
std::string symbol1 = "";
for (char ch: ucell.atoms[it].label)
{
if (std::isalpha(ch))
{
symbol1.push_back(ch);
}
}

double covalent_length
= (symbol1_covalent_radius + symbol2_covalent_radius)
/ ModuleBase::BOHR_TO_A;

for (int ia2 = 0; ia2 < ucell.atoms[it2].na; ia2++) {
for (int a = -1; a < 2; a++) {
for (int b = -1; b < 2; b++) {
for (int c = -1; c < 2; c++) {
if (it1 > it2) {
continue;
} else if (it1 == it2 && ia1 > ia2) {
continue;
} else if (it1 == it2 && ia1 == ia2 && a == 0
&& b == 0 && c == 0) {
continue;
}

double x2 = ucell.atoms[it2].taud[ia2].x + a;
double y2 = ucell.atoms[it2].taud[ia2].y + b;
double z2 = ucell.atoms[it2].taud[ia2].z + c;

double bond_length
= sqrt(pow((x2 - x1) * ucell.a1.x
+ (y2 - y1) * ucell.a2.x
+ (z2 - z1) * ucell.a3.x,
2)
+ pow((x2 - x1) * ucell.a1.y
+ (y2 - y1) * ucell.a2.y
+ (z2 - z1) * ucell.a3.y,
2)
+ pow((x2 - x1) * ucell.a1.z
+ (y2 - y1) * ucell.a2.z
+ (z2 - z1) * ucell.a3.z,
2))
* ucell.lat0;

if (bond_length < covalent_length * factor
|| bond_length
< covalent_length * warning_coef) {
errorlog.setf(std::ios_base::fixed,
std::ios_base::floatfield);
errorlog << std::setw(3) << ia1 + 1
<< "-th " << std::setw(3)
<< ucell.atoms[it1].label << ", ";
errorlog << std::setw(3) << ia2 + 1
<< "-th " << std::setw(3)
<< ucell.atoms[it2].label;
errorlog << " (cell:" << std::setw(2) << a
<< " " << std::setw(2) << b << " "
<< std::setw(2) << c << ")";
errorlog << ", distance= "
<< std::setprecision(3)
<< bond_length << " Bohr (";
errorlog
<< bond_length * ModuleBase::BOHR_TO_A
<< " Angstrom)" << std::endl;
if (ModuleBase::CovalentRadius.find(symbol1) != ModuleBase::CovalentRadius.end())
{
symbol_covalent_radiuss[it] = ModuleBase::CovalentRadius.at(symbol1);
}
else
{
std::stringstream mess;
mess << "Notice: symbol '" << symbol1 << "' is not an element symbol!!!! ";
mess << "set the covalent radius to be 0." << std::endl;
GlobalV::ofs_running << mess.str();
std::cout << mess.str();
}
}
std::vector<double> latvec (9);
latvec[0] = ucell.a1.x;
latvec[1] = ucell.a2.x;
latvec[2] = ucell.a3.x;
latvec[3] = ucell.a1.y;
latvec[4] = ucell.a2.y;
latvec[5] = ucell.a3.y;
latvec[6] = ucell.a1.z;
latvec[7] = ucell.a2.z;
latvec[8] = ucell.a3.z;
std::vector<double> A(27*3);
std::vector<std::string> cell(27);
std::vector<std::string> label(ntype);
for (int i = 0; i < 27; i++)
{
int a = (i / 9) % 3 - 1;
int b = (i / 3) % 3 - 1;
int c = i % 3 - 1;
A[3 * i] = a * latvec[0] + b * latvec[1] + c * latvec[2];
A[3 * i + 1] = a * latvec[3] + b * latvec[4] + c * latvec[5];
A[3 * i + 2] = a * latvec[6] + b * latvec[7] + c * latvec[8];
std::ostringstream tmp_oss;
tmp_oss << " (cell:" << std::setw(2) << a << " " << std::setw(2) << b << " " << std::setw(2) << c
<< "), distance= ";
cell[i] = tmp_oss.str();
}
for (int it = 0; it < ntype; it++)
{
std::ostringstream tmp_oss;
tmp_oss << std::setw(3) << ucell.atoms[it].label;
label[it] = tmp_oss.str();
}

if (bond_length
< covalent_length * factor) {
all_pass = false;
} else {
no_warning = false;
}
const double bohr_to_a = ModuleBase::BOHR_TO_A;
#pragma omp parallel
{
std::vector<double> delta_lat(3);
#pragma omp for schedule(dynamic)
for (int iat = 0; iat < ucell.nat; iat++)
{
const int it1 = ucell.iat2it[iat];
const int ia1 = ucell.iat2ia[iat];
const double symbol1_covalent_radius = symbol_covalent_radiuss[it1];
double x1 = ucell.atoms[it1].taud[ia1].x;
double y1 = ucell.atoms[it1].taud[ia1].y;
double z1 = ucell.atoms[it1].taud[ia1].z;
for (int it2 = it1; it2 < ntype; it2++)
{
double symbol2_covalent_radius = symbol_covalent_radiuss[it2];
double covalent_length = (symbol1_covalent_radius + symbol2_covalent_radius) / bohr_to_a;
const double max_error = covalent_length * max_factor_coef / ucell.lat0;
const double max_error_2 = max_error * max_error;
const double factor_error = covalent_length * factor;
for (int ia2 = ia1; ia2 < ucell.atoms[it2].na; ia2++)
{
const bool is_same_atom = (it1 == it2) && (ia1 == ia2);
double delta_x = ucell.atoms[it2].taud[ia2].x - x1;
double delta_y = ucell.atoms[it2].taud[ia2].y - y1;
double delta_z = ucell.atoms[it2].taud[ia2].z - z1;
delta_lat[0] = delta_x * latvec[0] + delta_y * latvec[1] + delta_z * latvec[2];
delta_lat[1] = delta_x * latvec[3] + delta_y * latvec[4] + delta_z * latvec[5];
delta_lat[2] = delta_x * latvec[6] + delta_y * latvec[7] + delta_z * latvec[8];
for (int i = 0; i < 27; i++)
{
if ((is_same_atom) && (i == 13))
continue;
const int offset = i * 3;
const double part1 = delta_lat[0] + A[offset];
const double part2 = delta_lat[1] + A[offset + 1];
const double part3 = delta_lat[2] + A[offset + 2];
const double bond_length = part1 * part1 + part2 * part2 + part3 * part3;
const bool flag = bond_length < max_error_2 ? true : false;
if (flag)
{
const double sqrt_bon = sqrt(bond_length) * lat0;
#pragma omp critical
{
no_warning = false;
all_pass = all_pass && (sqrt_bon < factor_error ? false : true);
errorlog << std::setw(3) << ia1 + 1 << "-th " << label[it1] << ", " << std::setw(3)
<< ia2 + 1 << "-th " << label[it2] << cell[i] << std::setprecision(3)
<< sqrt_bon << " Bohr (" << sqrt_bon * bohr_to_a << " Angstrom)\n";
}
} // c
} // b
} // a
} // ia2
} // it2
} // ia1
} // it1

if (!all_pass || !no_warning) {
}
}
} // ia2
} // it2
} // iat
}
ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru");
}
if (!all_pass || !no_warning)
{
std::stringstream mess;
mess << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
<< std::endl;
mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%"
<< std::endl;
mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
<< std::endl;
mess << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl;
mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl;
mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl;
mess << "!!! WARNING: Some atoms are too close!!!" << std::endl;
mess << "!!! Please check the nearest-neighbor list in log file."
<< std::endl;
mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
<< std::endl;
mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%"
<< std::endl;
mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
<< std::endl;
mess << "!!! Please check the nearest-neighbor list in log file." << std::endl;
mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl;
mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl;
mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl;

GlobalV::ofs_running << mess.str() << mess.str() << mess.str()
<< errorlog.str();
GlobalV::ofs_running << mess.str() << mess.str() << mess.str() << errorlog.str();
std::cout << mess.str() << mess.str() << mess.str() << std::endl;

if (!all_pass) {
if (!all_pass)
{
mess.clear();
mess.str("");
mess << "If this structure is what you want, you can set "
"'min_dist_coef'"
<< std::endl;
mess << "as a smaller value (the current value is " << factor
<< ") in INPUT file." << std::endl;
mess << "If this structure is what you want, you can set 'min_dist_coef'\n";
mess << "as a smaller value (the current value is " << factor << ") in INPUT file." << std::endl;
GlobalV::ofs_running << mess.str();
std::cout << mess.str();
ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!");
Expand Down
26 changes: 17 additions & 9 deletions source/module_cell/test/unitcell_test_readpp.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
#include "module_elecstate/read_pseudo.h"
#include <valarray>
#include <vector>
#include "string.h"
#ifdef __MPI
#include "mpi.h"
#endif
Expand Down Expand Up @@ -341,38 +342,45 @@ TEST_F(UcellDeathTest, CheckStructure) {
EXPECT_FALSE(ucell->atoms[0].ncpp.has_so);
EXPECT_FALSE(ucell->atoms[1].ncpp.has_so);
// trial 1

testing::internal::CaptureStdout();
double factor = 0.2;
ucell->set_iat2itia();
EXPECT_NO_THROW(Check_Atomic_Stru::check_atomic_stru(*ucell, factor));
output = testing::internal::GetCapturedStdout();
EXPECT_THAT(output,
testing::HasSubstr("WARNING: Some atoms are too close!!!"));
EXPECT_THAT(output,testing::HasSubstr("WARNING: Some atoms are too close!!!"));
// trial 2
testing::internal::CaptureStdout();
GlobalV::ofs_running.open("CheckStructure2.txt");
::testing::FLAGS_gtest_death_test_style = "threadsafe";
factor = 0.4;
EXPECT_EXIT(Check_Atomic_Stru::check_atomic_stru(*ucell, factor),
::testing::ExitedWithCode(1),
"");
output = testing::internal::GetCapturedStdout();
std::ifstream ifs("CheckStructure2.txt");
if (ifs.is_open())
{
std::string line;
while (std::getline(ifs, line)) {
output+=line;
}
}
EXPECT_THAT(output, testing::HasSubstr("The structure is unreasonable!"));
GlobalV::ofs_running.open("running.log");
// trial 3
ucell->atoms[0].label = "arbitrary";
testing::internal::CaptureStdout();
factor = 0.2;
EXPECT_NO_THROW(Check_Atomic_Stru::check_atomic_stru(*ucell, factor));
output = testing::internal::GetCapturedStdout();
EXPECT_THAT(
output,
testing::HasSubstr("Notice: symbol 'arbitrary' is not an element "
EXPECT_THAT(output,testing::HasSubstr("Notice: symbol 'arbitrary' is not an element "
"symbol!!!! set the covalent radius to be 0."));
// trial 4
ucell->atoms[0].label = "Fe1";
testing::internal::CaptureStdout();
factor = 0.2;
EXPECT_NO_THROW(Check_Atomic_Stru::check_atomic_stru(*ucell, factor));
output = testing::internal::GetCapturedStdout();
EXPECT_THAT(output,
testing::HasSubstr("WARNING: Some atoms are too close!!!"));
EXPECT_THAT(output,testing::HasSubstr("WARNING: Some atoms are too close!!!"));
}

TEST_F(UcellDeathTest, ReadPseudoWarning1) {
Expand Down
Loading