From a147dc28c6459c899358a65077053c8afbd734a9 Mon Sep 17 00:00:00 2001 From: tang070205 Date: Mon, 27 Jan 2025 22:57:46 +0800 Subject: [PATCH 01/20] add calculation of minimum atomic distance --- src/model/read_xyz.cu | 60 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/src/model/read_xyz.cu b/src/model/read_xyz.cu index 05e3f3f7c..e0f8c01a9 100644 --- a/src/model/read_xyz.cu +++ b/src/model/read_xyz.cu @@ -32,6 +32,7 @@ The class defining the simulation model. #include #include #include +#include const std::map MASS_TABLE{ {"H", 1.0080000000}, @@ -479,6 +480,64 @@ static std::vector get_atom_symbols(std::string& filename_potential return atom_symbols; } +static double calculate_distance( + const double* position_per_atom, int n1, int n2, int N, const Box& box) +{ + double dx = position_per_atom[n1] - position_per_atom[n2]; + double dy = position_per_atom[n1 + N] - position_per_atom[n2 + N]; + double dz = position_per_atom[n1 + 2 * N] - position_per_atom[n2 + 2 * N]; + + // 考虑周期性边界条件 + if (box.pbc_x) { + dx -= box.cpu_h[0] * std::round(box.cpu_h[9] * dx + box.cpu_h[10] * dy + box.cpu_h[11] * dz); + } + if (box.pbc_y) { + dy -= box.cpu_h[3] * std::round(box.cpu_h[12] * dx + box.cpu_h[13] * dy + box.cpu_h[14] * dz); + } + if (box.pbc_z) { + dz -= box.cpu_h[6] * std::round(box.cpu_h[15] * dx + box.cpu_h[16] * dy + box.cpu_h[17] * dz); + } + + return std::sqrt(dx * dx + dy * dy + dz * dz); +} + +static void check_min_distance(const Atom& atom, const Box& box) +{ + const int N = atom.number_of_atoms; + const double* position_per_atom = atom.cpu_position_per_atom.data(); + + double min_distance = 2.0; + int min_i = -1, min_j = -1; + + for (int i = 0; i < N; ++i) { + for (int j = i + 1; j < N; ++j) { + double distance = calculate_distance(position_per_atom, i, j, N, box); + if (distance < min_distance) { + min_distance = distance; + min_i = i; + min_j = j; + } + } + } + + if (min_distance < 2.0) { + printf("Error: Minimum distance (%f Å) between atoms %d (%s) and %d (%s) is less than 2 Å.\n", + min_distance, + min_i, + atom.cpu_atom_symbol[min_i].c_str(), + min_j, + atom.cpu_atom_symbol[min_j].c_str()); + PRINT_INPUT_ERROR("There are two atoms with a distance less than 2 Å."); + } else if (min_i != -1 && min_j != -1) { + printf("Minimum distance between atoms %d (%s) and %d (%s): %f Å\n", + min_i, + atom.cpu_atom_symbol[min_i].c_str(), + min_j, + atom.cpu_atom_symbol[min_j].c_str(), + min_distance); + } +} + void initialize_position( int& has_velocity_in_xyz, int& number_of_types, Box& box, std::vector& group, Atom& atom) { @@ -527,6 +586,7 @@ void initialize_position( } find_type_size(atom.number_of_atoms, number_of_types, atom.cpu_type, atom.cpu_type_size); + check_min_distance(atom, box); } void allocate_memory_gpu(std::vector& group, Atom& atom, GPU_Vector& thermo) From 89cf2751aa29699a4654ed72ef3dcb1b86f773b2 Mon Sep 17 00:00:00 2001 From: tang070205 Date: Mon, 27 Jan 2025 23:00:52 +0800 Subject: [PATCH 02/20] add calculation of minimum atomic distance --- src/model/read_xyz.cu | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/model/read_xyz.cu b/src/model/read_xyz.cu index e0f8c01a9..6348f3e61 100644 --- a/src/model/read_xyz.cu +++ b/src/model/read_xyz.cu @@ -520,14 +520,14 @@ static void check_min_distance(const Atom& atom, const Box& box) } } - if (min_distance < 2.0) { - printf("Error: Minimum distance (%f Å) between atoms %d (%s) and %d (%s) is less than 2 Å.\n", + if (min_distance < 1.0) { + printf("Error: Minimum distance (%f Å) between atoms %d (%s) and %d (%s) is less than 1 Å.\n", min_distance, min_i, atom.cpu_atom_symbol[min_i].c_str(), min_j, atom.cpu_atom_symbol[min_j].c_str()); - PRINT_INPUT_ERROR("There are two atoms with a distance less than 2 Å."); + PRINT_INPUT_ERROR("There are two atoms with a distance less than 1 Å."); } else if (min_i != -1 && min_j != -1) { printf("Minimum distance between atoms %d (%s) and %d (%s): %f Å\n", min_i, From d6766b0198fff2de1f4067e3cb087affa4f512e6 Mon Sep 17 00:00:00 2001 From: tang070205 Date: Wed, 29 Jan 2025 20:59:07 +0800 Subject: [PATCH 03/20] Add files via upload --- src/model/check_distance.cu | 143 +++++++++++++++++++++++++++++++++++ src/model/check_distance.cuh | 21 +++++ 2 files changed, 164 insertions(+) create mode 100644 src/model/check_distance.cu create mode 100644 src/model/check_distance.cuh diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu new file mode 100644 index 000000000..ce4e73bbc --- /dev/null +++ b/src/model/check_distance.cu @@ -0,0 +1,143 @@ +/* + Copyright 2017 Zheyong Fan and GPUMD development team + This file is part of GPUMD. + GPUMD is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + GPUMD is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with GPUMD. If not, see . +*/ + +/*----------------------------------------------------------------------------80 +Calculate the distance between any two atoms in the model.xyz file. +------------------------------------------------------------------------------*/ + +#include "check_distance.cuh" +#include "utilities/error.cuh" +#include "atom.cuh" +#include "box.cuh" +#include +#include +#include + +void calculate_min_atomic_distance(const Atom& atom, const Box& box) { + const int N = atom.number_of_atoms; + const double* pos = atom.cpu_position_per_atom.data(); + + double min_distance = 5.0; + int min_i = -1, min_j = -1; + + double Lx = sqrt(box.cpu_h[0]*box.cpu_h[0] + box.cpu_h[1]*box.cpu_h[1] + box.cpu_h[2]*box.cpu_h[2]); + double Ly = sqrt(box.cpu_h[3]*box.cpu_h[3] + box.cpu_h[4]*box.cpu_h[4] + box.cpu_h[5]*box.cpu_h[5]); + double Lz = sqrt(box.cpu_h[6]*box.cpu_h[6] + box.cpu_h[7]*box.cpu_h[7] + box.cpu_h[8]*box.cpu_h[8]); + + int nx = std::max(1, static_cast(ceil(Lx * 0.2))); + int ny = std::max(1, static_cast(ceil(Ly * 0.2))); + int nz = std::max(1, static_cast(ceil(Lz * 0.2))); + const double dx = Lx / nx, dy = Ly / ny, dz = Lz / nz; + + std::vector> wrapped_pos(N); + for (int i = 0; i < N; ++i) { + double x = pos[i]; + double y = pos[i + N]; + double z = pos[i + 2 * N]; + + double sx = box.cpu_h[9]*x + box.cpu_h[10]*y + box.cpu_h[11]*z; + double sy = box.cpu_h[12]*x + box.cpu_h[13]*y + box.cpu_h[14]*z; + double sz = box.cpu_h[15]*x + box.cpu_h[16]*y + box.cpu_h[17]*z; + + if (box.pbc_x) sx -= floor(sx); + if (box.pbc_y) sy -= floor(sy); + if (box.pbc_z) sz -= floor(sz); + + wrapped_pos[i][0] = sx*box.cpu_h[0] + sy*box.cpu_h[3] + sz*box.cpu_h[6]; + wrapped_pos[i][1] = sx*box.cpu_h[1] + sy*box.cpu_h[4] + sz*box.cpu_h[7]; + wrapped_pos[i][2] = sx*box.cpu_h[2] + sy*box.cpu_h[5] + sz*box.cpu_h[8]; + } + + std::vector> grid(nx * ny * nz); + + auto get_cell_index = [&](const std::array& p) { + int ix = static_cast((p[0] / Lx) * nx) % nx; + int iy = static_cast((p[1] / Ly) * ny) % ny; + int iz = static_cast((p[2] / Lz) * nz) % nz; + return ix + iy * nx + iz * nx * ny; + }; + + for (int i = 0; i < N; ++i) { + int cell = get_cell_index(wrapped_pos[i]); + grid[cell].push_back(i); + } + + for (int i = 0; i < N; ++i) { + const auto& pi = wrapped_pos[i]; + int cx = static_cast((pi[0]/Lx)*nx) % nx; + int cy = static_cast((pi[1]/Ly)*ny) % ny; + int cz = static_cast((pi[2]/Lz)*nz) % nz; + + for (int dx = -1; dx <= 1; ++dx) { + for (int dy = -1; dy <= 1; ++dy) { + for (int dz = -1; dz <= 1; ++dz) { + int nx_cell = (cx + dx + nx) % nx; + int ny_cell = (cy + dy + ny) % ny; + int nz_cell = (cz + dz + nz) % nz; + + if (!box.pbc_x && (cx + dx < 0 || cx + dx >= nx)) continue; + if (!box.pbc_y && (cy + dy < 0 || cy + dy >= ny)) continue; + if (!box.pbc_z && (cz + dz < 0 || cz + dz >= nz)) continue; + + int neighbor_cell = nx_cell + ny_cell * nx + nz_cell * nx * ny; + + for (int j : grid[neighbor_cell]) { + if (j <= i) continue; + + const auto& pj = wrapped_pos[j]; + double delta[3] = {pi[0]-pj[0], pi[1]-pj[1], pi[2]-pj[2]}; + + if (fabs(delta[0]) > 2.0 || fabs(delta[1]) > 2.0 || fabs(delta[2]) > 2.0) continue; + + if (box.pbc_x) {delta[0] -= box.cpu_h[0] * std::round(box.cpu_h[9] * delta[0] + box.cpu_h[10] * delta[1] + box.cpu_h[11] * delta[2]); + } + if (box.pbc_y) {delta[1] -= box.cpu_h[3] * std::round(box.cpu_h[12] * delta[0] + box.cpu_h[13] * delta[1] + box.cpu_h[14] * delta[2]); + } + if (box.pbc_z) {delta[2] -= box.cpu_h[6] * std::round(box.cpu_h[15] * delta[0] + box.cpu_h[16] * delta[1] + box.cpu_h[17] * delta[2]); + } + + double dist_sq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2]; + if (dist_sq >= 4.0) continue; + + double dist = sqrt(dist_sq); + if (dist < min_distance) { + min_distance = dist; + min_i = i; + min_j = j; + } + } + } + } + } + } + + if (min_distance < 1.0) { + printf("Error: Minimum distance (%f Å) between atoms %d (%s) and %d (%s) is less than 1 Å.\n", + min_distance, + min_i, + atom.cpu_atom_symbol[min_i].c_str(), + min_j, + atom.cpu_atom_symbol[min_j].c_str()); + PRINT_INPUT_ERROR("There are two atoms with a distance less than 1 Å."); + } else if (min_i != -1 && min_j != -1) { + printf("Minimum distance between atoms %d (%s) and %d (%s): %f Å\n", + min_i, + atom.cpu_atom_symbol[min_i].c_str(), + min_j, + atom.cpu_atom_symbol[min_j].c_str(), + min_distance); + } +} + diff --git a/src/model/check_distance.cuh b/src/model/check_distance.cuh new file mode 100644 index 000000000..903602b76 --- /dev/null +++ b/src/model/check_distance.cuh @@ -0,0 +1,21 @@ +/* + Copyright 2017 Zheyong Fan and GPUMD development team + This file is part of GPUMD. + GPUMD is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + GPUMD is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with GPUMD. If not, see . +*/ + +#pragma once + +class Box; +class Atom; + +void calculate_min_atomic_distance(const Atom& atom, const Box& box); \ No newline at end of file From e58e9935b3bf5b51279b08b4fae9461419341152 Mon Sep 17 00:00:00 2001 From: tang070205 Date: Wed, 29 Jan 2025 21:00:36 +0800 Subject: [PATCH 04/20] Update run.cu --- src/main_gpumd/run.cu | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/main_gpumd/run.cu b/src/main_gpumd/run.cu index dc714b6f9..b99101cb1 100644 --- a/src/main_gpumd/run.cu +++ b/src/main_gpumd/run.cu @@ -29,6 +29,7 @@ Run simulation according to the inputs in the run.in file. #include "minimize/minimize.cuh" #include "model/box.cuh" #include "model/read_xyz.cuh" +#include "model/check_distance.cuh" #include "phonon/hessian.cuh" #include "replicate.cuh" #include "run.cuh" @@ -115,6 +116,8 @@ Run::Run() initialize_position(has_velocity_in_xyz, number_of_types, box, group, atom); + calculate_min_atomic_distance(atom, box); + allocate_memory_gpu(group, atom, thermo); velocity.initialize( From e669ff73e52cb768ba438939a117038d08a148ed Mon Sep 17 00:00:00 2001 From: tang070205 Date: Wed, 29 Jan 2025 21:02:45 +0800 Subject: [PATCH 05/20] Update read_xyz.cu --- src/model/read_xyz.cu | 60 ------------------------------------------- 1 file changed, 60 deletions(-) diff --git a/src/model/read_xyz.cu b/src/model/read_xyz.cu index 6348f3e61..05e3f3f7c 100644 --- a/src/model/read_xyz.cu +++ b/src/model/read_xyz.cu @@ -32,7 +32,6 @@ The class defining the simulation model. #include #include #include -#include const std::map MASS_TABLE{ {"H", 1.0080000000}, @@ -480,64 +479,6 @@ static std::vector get_atom_symbols(std::string& filename_potential return atom_symbols; } -static double calculate_distance( - const double* position_per_atom, int n1, int n2, int N, const Box& box) -{ - double dx = position_per_atom[n1] - position_per_atom[n2]; - double dy = position_per_atom[n1 + N] - position_per_atom[n2 + N]; - double dz = position_per_atom[n1 + 2 * N] - position_per_atom[n2 + 2 * N]; - - // 考虑周期性边界条件 - if (box.pbc_x) { - dx -= box.cpu_h[0] * std::round(box.cpu_h[9] * dx + box.cpu_h[10] * dy + box.cpu_h[11] * dz); - } - if (box.pbc_y) { - dy -= box.cpu_h[3] * std::round(box.cpu_h[12] * dx + box.cpu_h[13] * dy + box.cpu_h[14] * dz); - } - if (box.pbc_z) { - dz -= box.cpu_h[6] * std::round(box.cpu_h[15] * dx + box.cpu_h[16] * dy + box.cpu_h[17] * dz); - } - - return std::sqrt(dx * dx + dy * dy + dz * dz); -} - -static void check_min_distance(const Atom& atom, const Box& box) -{ - const int N = atom.number_of_atoms; - const double* position_per_atom = atom.cpu_position_per_atom.data(); - - double min_distance = 2.0; - int min_i = -1, min_j = -1; - - for (int i = 0; i < N; ++i) { - for (int j = i + 1; j < N; ++j) { - double distance = calculate_distance(position_per_atom, i, j, N, box); - if (distance < min_distance) { - min_distance = distance; - min_i = i; - min_j = j; - } - } - } - - if (min_distance < 1.0) { - printf("Error: Minimum distance (%f Å) between atoms %d (%s) and %d (%s) is less than 1 Å.\n", - min_distance, - min_i, - atom.cpu_atom_symbol[min_i].c_str(), - min_j, - atom.cpu_atom_symbol[min_j].c_str()); - PRINT_INPUT_ERROR("There are two atoms with a distance less than 1 Å."); - } else if (min_i != -1 && min_j != -1) { - printf("Minimum distance between atoms %d (%s) and %d (%s): %f Å\n", - min_i, - atom.cpu_atom_symbol[min_i].c_str(), - min_j, - atom.cpu_atom_symbol[min_j].c_str(), - min_distance); - } -} - void initialize_position( int& has_velocity_in_xyz, int& number_of_types, Box& box, std::vector& group, Atom& atom) { @@ -586,7 +527,6 @@ void initialize_position( } find_type_size(atom.number_of_atoms, number_of_types, atom.cpu_type, atom.cpu_type_size); - check_min_distance(atom, box); } void allocate_memory_gpu(std::vector& group, Atom& atom, GPU_Vector& thermo) From ac4862f3656a76bbfa5e9a63e30df4e1e198e9ee Mon Sep 17 00:00:00 2001 From: tang070205 Date: Wed, 29 Jan 2025 22:50:15 +0800 Subject: [PATCH 06/20] Update check_distance.cu --- src/model/check_distance.cu | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index ce4e73bbc..162a3c49b 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -32,9 +32,9 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) { double min_distance = 5.0; int min_i = -1, min_j = -1; - double Lx = sqrt(box.cpu_h[0]*box.cpu_h[0] + box.cpu_h[1]*box.cpu_h[1] + box.cpu_h[2]*box.cpu_h[2]); - double Ly = sqrt(box.cpu_h[3]*box.cpu_h[3] + box.cpu_h[4]*box.cpu_h[4] + box.cpu_h[5]*box.cpu_h[5]); - double Lz = sqrt(box.cpu_h[6]*box.cpu_h[6] + box.cpu_h[7]*box.cpu_h[7] + box.cpu_h[8]*box.cpu_h[8]); + double Lx = sqrt(box.cpu_h[0]*box.cpu_h[0] + box.cpu_h[3]*box.cpu_h[3] + box.cpu_h[6]*box.cpu_h[6]); + double Ly = sqrt(box.cpu_h[1]*box.cpu_h[1] + box.cpu_h[4]*box.cpu_h[4] + box.cpu_h[7]*box.cpu_h[7]); + double Lz = sqrt(box.cpu_h[2]*box.cpu_h[2] + box.cpu_h[5]*box.cpu_h[5] + box.cpu_h[8]*box.cpu_h[8]); int nx = std::max(1, static_cast(ceil(Lx * 0.2))); int ny = std::max(1, static_cast(ceil(Ly * 0.2))); From eec3f3d35350e839e82424ee5943f9dbfaf710a7 Mon Sep 17 00:00:00 2001 From: tang070205 Date: Wed, 29 Jan 2025 23:44:48 +0800 Subject: [PATCH 07/20] Update check_distance.cu --- src/model/check_distance.cu | 249 ++++++++++++++++++++---------------- 1 file changed, 136 insertions(+), 113 deletions(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index 162a3c49b..17c8aafe4 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -17,127 +17,150 @@ Calculate the distance between any two atoms in the model.xyz file. ------------------------------------------------------------------------------*/ -#include "check_distance.cuh" -#include "utilities/error.cuh" #include "atom.cuh" #include "box.cuh" -#include +#include "check_distance.cuh" +#include "utilities/error.cuh" #include +#include #include -void calculate_min_atomic_distance(const Atom& atom, const Box& box) { - const int N = atom.number_of_atoms; - const double* pos = atom.cpu_position_per_atom.data(); - - double min_distance = 5.0; - int min_i = -1, min_j = -1; - - double Lx = sqrt(box.cpu_h[0]*box.cpu_h[0] + box.cpu_h[3]*box.cpu_h[3] + box.cpu_h[6]*box.cpu_h[6]); - double Ly = sqrt(box.cpu_h[1]*box.cpu_h[1] + box.cpu_h[4]*box.cpu_h[4] + box.cpu_h[7]*box.cpu_h[7]); - double Lz = sqrt(box.cpu_h[2]*box.cpu_h[2] + box.cpu_h[5]*box.cpu_h[5] + box.cpu_h[8]*box.cpu_h[8]); - - int nx = std::max(1, static_cast(ceil(Lx * 0.2))); - int ny = std::max(1, static_cast(ceil(Ly * 0.2))); - int nz = std::max(1, static_cast(ceil(Lz * 0.2))); - const double dx = Lx / nx, dy = Ly / ny, dz = Lz / nz; - - std::vector> wrapped_pos(N); - for (int i = 0; i < N; ++i) { - double x = pos[i]; - double y = pos[i + N]; - double z = pos[i + 2 * N]; - - double sx = box.cpu_h[9]*x + box.cpu_h[10]*y + box.cpu_h[11]*z; - double sy = box.cpu_h[12]*x + box.cpu_h[13]*y + box.cpu_h[14]*z; - double sz = box.cpu_h[15]*x + box.cpu_h[16]*y + box.cpu_h[17]*z; - - if (box.pbc_x) sx -= floor(sx); - if (box.pbc_y) sy -= floor(sy); - if (box.pbc_z) sz -= floor(sz); - - wrapped_pos[i][0] = sx*box.cpu_h[0] + sy*box.cpu_h[3] + sz*box.cpu_h[6]; - wrapped_pos[i][1] = sx*box.cpu_h[1] + sy*box.cpu_h[4] + sz*box.cpu_h[7]; - wrapped_pos[i][2] = sx*box.cpu_h[2] + sy*box.cpu_h[5] + sz*box.cpu_h[8]; - } - - std::vector> grid(nx * ny * nz); - - auto get_cell_index = [&](const std::array& p) { - int ix = static_cast((p[0] / Lx) * nx) % nx; - int iy = static_cast((p[1] / Ly) * ny) % ny; - int iz = static_cast((p[2] / Lz) * nz) % nz; - return ix + iy * nx + iz * nx * ny; - }; +void calculate_min_atomic_distance(const Atom& atom, const Box& box) +{ + const int N = atom.number_of_atoms; + const double* pos = atom.cpu_position_per_atom.data(); + + double min_distance = 5.0; + int min_i = -1, min_j = -1; + + double Lx = + sqrt(box.cpu_h[0] * box.cpu_h[0] + box.cpu_h[1] * box.cpu_h[1] + box.cpu_h[2] * box.cpu_h[2]); + double Ly = + sqrt(box.cpu_h[3] * box.cpu_h[3] + box.cpu_h[4] * box.cpu_h[4] + box.cpu_h[5] * box.cpu_h[5]); + double Lz = + sqrt(box.cpu_h[6] * box.cpu_h[6] + box.cpu_h[7] * box.cpu_h[7] + box.cpu_h[8] * box.cpu_h[8]); + + int nx = std::max(1, static_cast(ceil(Lx * 0.2))); + int ny = std::max(1, static_cast(ceil(Ly * 0.2))); + int nz = std::max(1, static_cast(ceil(Lz * 0.2))); + const double dx = Lx / nx, dy = Ly / ny, dz = Lz / nz; + + std::vector> wrapped_pos(N); + for (int i = 0; i < N; ++i) { + double x = pos[i]; + double y = pos[i + N]; + double z = pos[i + 2 * N]; + + double sx = box.cpu_h[9] * x + box.cpu_h[10] * y + box.cpu_h[11] * z; + double sy = box.cpu_h[12] * x + box.cpu_h[13] * y + box.cpu_h[14] * z; + double sz = box.cpu_h[15] * x + box.cpu_h[16] * y + box.cpu_h[17] * z; + + if (box.pbc_x) + sx -= floor(sx); + if (box.pbc_y) + sy -= floor(sy); + if (box.pbc_z) + sz -= floor(sz); + + wrapped_pos[i][0] = sx * box.cpu_h[0] + sy * box.cpu_h[3] + sz * box.cpu_h[6]; + wrapped_pos[i][1] = sx * box.cpu_h[1] + sy * box.cpu_h[4] + sz * box.cpu_h[7]; + wrapped_pos[i][2] = sx * box.cpu_h[2] + sy * box.cpu_h[5] + sz * box.cpu_h[8]; + } + + std::vector> grid(nx * ny * nz); + + auto get_cell_index = [&](const std::array& p) { + int ix = static_cast((p[0] / Lx) * nx) % nx; + int iy = static_cast((p[1] / Ly) * ny) % ny; + int iz = static_cast((p[2] / Lz) * nz) % nz; + return ix + iy * nx + iz * nx * ny; + }; + + for (int i = 0; i < N; ++i) { + int cell = get_cell_index(wrapped_pos[i]); + grid[cell].push_back(i); + } + + for (int i = 0; i < N; ++i) { + const auto& pi = wrapped_pos[i]; + int cx = static_cast((pi[0] / Lx) * nx) % nx; + int cy = static_cast((pi[1] / Ly) * ny) % ny; + int cz = static_cast((pi[2] / Lz) * nz) % nz; + + for (int dx = -1; dx <= 1; ++dx) { + for (int dy = -1; dy <= 1; ++dy) { + for (int dz = -1; dz <= 1; ++dz) { + int nx_cell = (cx + dx + nx) % nx; + int ny_cell = (cy + dy + ny) % ny; + int nz_cell = (cz + dz + nz) % nz; + + if (!box.pbc_x && (cx + dx < 0 || cx + dx >= nx)) + continue; + if (!box.pbc_y && (cy + dy < 0 || cy + dy >= ny)) + continue; + if (!box.pbc_z && (cz + dz < 0 || cz + dz >= nz)) + continue; + + int neighbor_cell = nx_cell + ny_cell * nx + nz_cell * nx * ny; + + for (int j : grid[neighbor_cell]) { + if (j <= i) + continue; + + const auto& pj = wrapped_pos[j]; + double delta[3] = {pi[0] - pj[0], pi[1] - pj[1], pi[2] - pj[2]}; + + if (fabs(delta[0]) > 2.0 || fabs(delta[1]) > 2.0 || fabs(delta[2]) > 2.0) + continue; + + if (box.pbc_x) { + delta[0] -= box.cpu_h[0] * std::round( + box.cpu_h[9] * delta[0] + box.cpu_h[10] * delta[1] + + box.cpu_h[11] * delta[2]); + } + if (box.pbc_y) { + delta[1] -= box.cpu_h[3] * std::round( + box.cpu_h[12] * delta[0] + box.cpu_h[13] * delta[1] + + box.cpu_h[14] * delta[2]); + } + if (box.pbc_z) { + delta[2] -= box.cpu_h[6] * std::round( + box.cpu_h[15] * delta[0] + box.cpu_h[16] * delta[1] + + box.cpu_h[17] * delta[2]); + } - for (int i = 0; i < N; ++i) { - int cell = get_cell_index(wrapped_pos[i]); - grid[cell].push_back(i); - } + double dist_sq = delta[0] * delta[0] + delta[1] * delta[1] + delta[2] * delta[2]; + if (dist_sq >= 4.0) + continue; - for (int i = 0; i < N; ++i) { - const auto& pi = wrapped_pos[i]; - int cx = static_cast((pi[0]/Lx)*nx) % nx; - int cy = static_cast((pi[1]/Ly)*ny) % ny; - int cz = static_cast((pi[2]/Lz)*nz) % nz; - - for (int dx = -1; dx <= 1; ++dx) { - for (int dy = -1; dy <= 1; ++dy) { - for (int dz = -1; dz <= 1; ++dz) { - int nx_cell = (cx + dx + nx) % nx; - int ny_cell = (cy + dy + ny) % ny; - int nz_cell = (cz + dz + nz) % nz; - - if (!box.pbc_x && (cx + dx < 0 || cx + dx >= nx)) continue; - if (!box.pbc_y && (cy + dy < 0 || cy + dy >= ny)) continue; - if (!box.pbc_z && (cz + dz < 0 || cz + dz >= nz)) continue; - - int neighbor_cell = nx_cell + ny_cell * nx + nz_cell * nx * ny; - - for (int j : grid[neighbor_cell]) { - if (j <= i) continue; - - const auto& pj = wrapped_pos[j]; - double delta[3] = {pi[0]-pj[0], pi[1]-pj[1], pi[2]-pj[2]}; - - if (fabs(delta[0]) > 2.0 || fabs(delta[1]) > 2.0 || fabs(delta[2]) > 2.0) continue; - - if (box.pbc_x) {delta[0] -= box.cpu_h[0] * std::round(box.cpu_h[9] * delta[0] + box.cpu_h[10] * delta[1] + box.cpu_h[11] * delta[2]); - } - if (box.pbc_y) {delta[1] -= box.cpu_h[3] * std::round(box.cpu_h[12] * delta[0] + box.cpu_h[13] * delta[1] + box.cpu_h[14] * delta[2]); - } - if (box.pbc_z) {delta[2] -= box.cpu_h[6] * std::round(box.cpu_h[15] * delta[0] + box.cpu_h[16] * delta[1] + box.cpu_h[17] * delta[2]); - } - - double dist_sq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2]; - if (dist_sq >= 4.0) continue; - - double dist = sqrt(dist_sq); - if (dist < min_distance) { - min_distance = dist; - min_i = i; - min_j = j; - } - } - } + double dist = sqrt(dist_sq); + if (dist < min_distance) { + min_distance = dist; + min_i = i; + min_j = j; } + } } + } } - - if (min_distance < 1.0) { - printf("Error: Minimum distance (%f Å) between atoms %d (%s) and %d (%s) is less than 1 Å.\n", - min_distance, - min_i, - atom.cpu_atom_symbol[min_i].c_str(), - min_j, - atom.cpu_atom_symbol[min_j].c_str()); - PRINT_INPUT_ERROR("There are two atoms with a distance less than 1 Å."); - } else if (min_i != -1 && min_j != -1) { - printf("Minimum distance between atoms %d (%s) and %d (%s): %f Å\n", - min_i, - atom.cpu_atom_symbol[min_i].c_str(), - min_j, - atom.cpu_atom_symbol[min_j].c_str(), - min_distance); - } + } + + if (min_distance < 1.0) { + printf( + "Error: Minimum distance (%f Å) between atoms %d (%s) and %d (%s) is less than 1 Å.\n", + min_distance, + min_i, + atom.cpu_atom_symbol[min_i].c_str(), + min_j, + atom.cpu_atom_symbol[min_j].c_str()); + PRINT_INPUT_ERROR("There are two atoms with a distance less than 1 Å."); + } else if (min_i != -1 && min_j != -1) { + printf( + "Minimum distance between atoms %d (%s) and %d (%s): %f Å\n", + min_i, + atom.cpu_atom_symbol[min_i].c_str(), + min_j, + atom.cpu_atom_symbol[min_j].c_str(), + min_distance); + } } - From d6e6c748c5abde5c5d8de0d45de81f4f8350595d Mon Sep 17 00:00:00 2001 From: tang070205 Date: Thu, 30 Jan 2025 21:11:48 +0800 Subject: [PATCH 08/20] Update check_distance.cu --- src/model/check_distance.cu | 234 +++++++++++++++++++----------------- 1 file changed, 122 insertions(+), 112 deletions(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index 17c8aafe4..9cc440304 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -21,9 +21,51 @@ Calculate the distance between any two atoms in the model.xyz file. #include "box.cuh" #include "check_distance.cuh" #include "utilities/error.cuh" -#include #include -#include + +void applyMicOne(double& x12) +{ + if (x12 < -0.5) { + x12 += 1.0; + } else if (x12 > +0.5) { + x12 -= 1.0; + } +} + +void applyMic(const Box& box, double& x12, double& y12, double& z12) +{ + + double sx = box.cpu_h[9] * x12 + box.cpu_h[10] * y12 + box.cpu_h[11] * z12; + double sy = box.cpu_h[12] * x12 + box.cpu_h[13] * y12 + box.cpu_h[14] * z12; + double sz = box.cpu_h[15] * x12 + box.cpu_h[16] * y12 + box.cpu_h[17] * z12; + + if (box.pbc_x) + applyMicOne(sx); + if (box.pbc_y) + applyMicOne(sy); + if (box.pbc_z) + applyMicOne(sz); + + x12 = box.cpu_h[0] * sx + box.cpu_h[3] * sy + box.cpu_h[6] * sz; + y12 = box.cpu_h[1] * sx + box.cpu_h[4] * sy + box.cpu_h[7] * sz; + z12 = box.cpu_h[2] * sx + box.cpu_h[5] * sy + box.cpu_h[8] * sz; +} + +void findCell(const Box& box, const double* thickness, const double* r, const int* numCells, int* cell) +{ + double s[3]; + s[0] = box.cpu_h[9] * r[0] + box.cpu_h[10] * r[1] + box.cpu_h[11] * r[2]; + s[1] = box.cpu_h[12] * r[0] + box.cpu_h[13] * r[1] + box.cpu_h[14] * r[2]; + s[2] = box.cpu_h[15] * r[0] + box.cpu_h[16] * r[1] + box.cpu_h[17] * r[2]; + for (int d = 0; d < 3; ++d) { + cell[d] = floor(s[d] * thickness[d] * 0.2); + if (cell[d] < 0) + cell[d] += numCells[d]; + if (cell[d] >= numCells[d]) + cell[d] -= numCells[d]; + } + cell[3] = cell[0] + numCells[0] * (cell[1] + numCells[1] * cell[2]); +} void calculate_min_atomic_distance(const Atom& atom, const Box& box) { @@ -31,113 +73,81 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) const double* pos = atom.cpu_position_per_atom.data(); double min_distance = 5.0; - int min_i = -1, min_j = -1; - - double Lx = - sqrt(box.cpu_h[0] * box.cpu_h[0] + box.cpu_h[1] * box.cpu_h[1] + box.cpu_h[2] * box.cpu_h[2]); - double Ly = - sqrt(box.cpu_h[3] * box.cpu_h[3] + box.cpu_h[4] * box.cpu_h[4] + box.cpu_h[5] * box.cpu_h[5]); - double Lz = - sqrt(box.cpu_h[6] * box.cpu_h[6] + box.cpu_h[7] * box.cpu_h[7] + box.cpu_h[8] * box.cpu_h[8]); - - int nx = std::max(1, static_cast(ceil(Lx * 0.2))); - int ny = std::max(1, static_cast(ceil(Ly * 0.2))); - int nz = std::max(1, static_cast(ceil(Lz * 0.2))); - const double dx = Lx / nx, dy = Ly / ny, dz = Lz / nz; - - std::vector> wrapped_pos(N); - for (int i = 0; i < N; ++i) { - double x = pos[i]; - double y = pos[i + N]; - double z = pos[i + 2 * N]; - - double sx = box.cpu_h[9] * x + box.cpu_h[10] * y + box.cpu_h[11] * z; - double sy = box.cpu_h[12] * x + box.cpu_h[13] * y + box.cpu_h[14] * z; - double sz = box.cpu_h[15] * x + box.cpu_h[16] * y + box.cpu_h[17] * z; - - if (box.pbc_x) - sx -= floor(sx); - if (box.pbc_y) - sy -= floor(sy); - if (box.pbc_z) - sz -= floor(sz); - - wrapped_pos[i][0] = sx * box.cpu_h[0] + sy * box.cpu_h[3] + sz * box.cpu_h[6]; - wrapped_pos[i][1] = sx * box.cpu_h[1] + sy * box.cpu_h[4] + sz * box.cpu_h[7]; - wrapped_pos[i][2] = sx * box.cpu_h[2] + sy * box.cpu_h[5] + sz * box.cpu_h[8]; - } + int min_n1 = -1, min_n2 = -1; - std::vector> grid(nx * ny * nz); + double thickness[3]; + thickness[0] = sqrt(box.cpu_h[0] * box.cpu_h[0] + box.cpu_h[3] * box.cpu_h[3] + box.cpu_h[6] * box.cpu_h[6]); + thickness[1] = sqrt(box.cpu_h[1] * box.cpu_h[1] + box.cpu_h[4] * box.cpu_h[4] + box.cpu_h[7] * box.cpu_h[7]); + thickness[2] = sqrt(box.cpu_h[2] * box.cpu_h[2] + box.cpu_h[5] * box.cpu_h[5] + box.cpu_h[8] * box.cpu_h[8]); - auto get_cell_index = [&](const std::array& p) { - int ix = static_cast((p[0] / Lx) * nx) % nx; - int iy = static_cast((p[1] / Ly) * ny) % ny; - int iz = static_cast((p[2] / Lz) * nz) % nz; - return ix + iy * nx + iz * nx * ny; - }; + int numCells[4]; + numCells[0] = std::max(1, static_cast(ceil(thickness[0] * 0.2))); + numCells[1] = std::max(1, static_cast(ceil(thickness[1] * 0.2))); + numCells[2] = std::max(1, static_cast(ceil(thickness[2] * 0.2))); + numCells[3] = numCells[0] * numCells[1] * numCells[2]; - for (int i = 0; i < N; ++i) { - int cell = get_cell_index(wrapped_pos[i]); - grid[cell].push_back(i); + int cell[4]; + + std::vector cellCount(numCells[3], 0); + std::vector cellCountSum(numCells[3], 0); + + for (int n = 0; n < N; ++n) { + const double r[3] = {pos[n], pos[n + N], pos[n + 2 * N]}; + findCell(box, thickness, r, numCells, cell); + ++cellCount[cell[3]]; } - for (int i = 0; i < N; ++i) { - const auto& pi = wrapped_pos[i]; - int cx = static_cast((pi[0] / Lx) * nx) % nx; - int cy = static_cast((pi[1] / Ly) * ny) % ny; - int cz = static_cast((pi[2] / Lz) * nz) % nz; - - for (int dx = -1; dx <= 1; ++dx) { - for (int dy = -1; dy <= 1; ++dy) { - for (int dz = -1; dz <= 1; ++dz) { - int nx_cell = (cx + dx + nx) % nx; - int ny_cell = (cy + dy + ny) % ny; - int nz_cell = (cz + dz + nz) % nz; - - if (!box.pbc_x && (cx + dx < 0 || cx + dx >= nx)) - continue; - if (!box.pbc_y && (cy + dy < 0 || cy + dy >= ny)) - continue; - if (!box.pbc_z && (cz + dz < 0 || cz + dz >= nz)) - continue; - - int neighbor_cell = nx_cell + ny_cell * nx + nz_cell * nx * ny; - - for (int j : grid[neighbor_cell]) { - if (j <= i) - continue; - - const auto& pj = wrapped_pos[j]; - double delta[3] = {pi[0] - pj[0], pi[1] - pj[1], pi[2] - pj[2]}; - - if (fabs(delta[0]) > 2.0 || fabs(delta[1]) > 2.0 || fabs(delta[2]) > 2.0) - continue; - - if (box.pbc_x) { - delta[0] -= box.cpu_h[0] * std::round( - box.cpu_h[9] * delta[0] + box.cpu_h[10] * delta[1] + - box.cpu_h[11] * delta[2]); - } - if (box.pbc_y) { - delta[1] -= box.cpu_h[3] * std::round( - box.cpu_h[12] * delta[0] + box.cpu_h[13] * delta[1] + - box.cpu_h[14] * delta[2]); - } - if (box.pbc_z) { - delta[2] -= box.cpu_h[6] * std::round( - box.cpu_h[15] * delta[0] + box.cpu_h[16] * delta[1] + - box.cpu_h[17] * delta[2]); - } + for (int i = 1; i < numCells[3]; ++i) { + cellCountSum[i] = cellCountSum[i - 1] + cellCount[i - 1]; + } - double dist_sq = delta[0] * delta[0] + delta[1] * delta[1] + delta[2] * delta[2]; - if (dist_sq >= 4.0) - continue; + std::fill(cellCount.begin(), cellCount.end(), 0); + std::vector cellContents(N, 0); + + for (int n = 0; n < N; ++n) { + const double r[3] = {pos[n], pos[n + N], pos[n + 2 * N]}; + findCell(box, thickness, r, numCells, cell); + cellContents[cellCountSum[cell[3]] + cellCount[cell[3]]] = n; + ++cellCount[cell[3]]; + } - double dist = sqrt(dist_sq); - if (dist < min_distance) { - min_distance = dist; - min_i = i; - min_j = j; + for (int n1 = 0; n1 < N; ++n1) { + const double r1[3] = {pos[n1], pos[n1 + N], pos[n1 + 2 * N]}; + findCell(box, thickness, r1, numCells, cell); + for (int k = -1; k <= 1; ++k) { + for (int j = -1; j <= 1; ++j) { + for (int i = -1; i <= 1; ++i) { + int neighborCell = cell[3] + (k * numCells[1] + j) * numCells[0] + i; + if (cell[0] + i < 0) + neighborCell += numCells[0]; + if (cell[0] + i >= numCells[0]) + neighborCell -= numCells[0]; + if (cell[1] + j < 0) + neighborCell += numCells[1] * numCells[0]; + if (cell[1] + j >= numCells[1]) + neighborCell -= numCells[1] * numCells[0]; + if (cell[2] + k < 0) + neighborCell += numCells[3]; + if (cell[2] + k >= numCells[2]) + neighborCell -= numCells[3]; + for (int m = 0; m < cellCount[neighborCell]; ++m) { + const int n2 = cellContents[cellCountSum[neighborCell] + m]; + if (n1 < n2) { + double x12 = pos[n2] - r1[0]; + double y12 = pos[n2 + N] - r1[1]; + double z12 = pos[n2 + 2 * N] - r1[2]; + applyMic(box, x12, y12, z12); + if (fabs(x12) > 2.0 || fabs(y12) > 2.0 || fabs(z12) > 2.0) continue; + const double d2 = x12 * x12 + y12 * y12 + z12 * z12; + if (d2 >= 4.0) + continue; + + double dist = sqrt(d2); + if (dist < min_distance) { + min_distance = dist; + min_n1 = n1; + min_n2 = n2; + } } } } @@ -145,22 +155,22 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) } } - if (min_distance < 1.0) { + if (min_distance < 2.0) { printf( "Error: Minimum distance (%f Å) between atoms %d (%s) and %d (%s) is less than 1 Å.\n", min_distance, - min_i, - atom.cpu_atom_symbol[min_i].c_str(), - min_j, - atom.cpu_atom_symbol[min_j].c_str()); + min_n1, + atom.cpu_atom_symbol[min_n1].c_str(), + min_n2, + atom.cpu_atom_symbol[min_n2].c_str()); PRINT_INPUT_ERROR("There are two atoms with a distance less than 1 Å."); - } else if (min_i != -1 && min_j != -1) { + } else if (min_n1 != -1 && min_n2 != -1) { printf( "Minimum distance between atoms %d (%s) and %d (%s): %f Å\n", - min_i, - atom.cpu_atom_symbol[min_i].c_str(), - min_j, - atom.cpu_atom_symbol[min_j].c_str(), + min_n1, + atom.cpu_atom_symbol[min_n1].c_str(), + min_n2, + atom.cpu_atom_symbol[min_n2].c_str(), min_distance); } } From f4bded97f17ba9ccd2d7612f567fd3f9aabc80c2 Mon Sep 17 00:00:00 2001 From: tang070205 Date: Thu, 30 Jan 2025 21:13:41 +0800 Subject: [PATCH 09/20] Update check_distance.cu --- src/model/check_distance.cu | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index 9cc440304..98e63e81f 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -51,7 +51,8 @@ void applyMic(const Box& box, double& x12, double& y12, double& z12) z12 = box.cpu_h[2] * sx + box.cpu_h[5] * sy + box.cpu_h[8] * sz; } -void findCell(const Box& box, const double* thickness, const double* r, const int* numCells, int* cell) +void findCell( + const Box& box, const double* thickness, const double* r, const int* numCells, int* cell) { double s[3]; s[0] = box.cpu_h[9] * r[0] + box.cpu_h[10] * r[1] + box.cpu_h[11] * r[2]; @@ -76,9 +77,12 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) int min_n1 = -1, min_n2 = -1; double thickness[3]; - thickness[0] = sqrt(box.cpu_h[0] * box.cpu_h[0] + box.cpu_h[3] * box.cpu_h[3] + box.cpu_h[6] * box.cpu_h[6]); - thickness[1] = sqrt(box.cpu_h[1] * box.cpu_h[1] + box.cpu_h[4] * box.cpu_h[4] + box.cpu_h[7] * box.cpu_h[7]); - thickness[2] = sqrt(box.cpu_h[2] * box.cpu_h[2] + box.cpu_h[5] * box.cpu_h[5] + box.cpu_h[8] * box.cpu_h[8]); + thickness[0] = + sqrt(box.cpu_h[0] * box.cpu_h[0] + box.cpu_h[3] * box.cpu_h[3] + box.cpu_h[6] * box.cpu_h[6]); + thickness[1] = + sqrt(box.cpu_h[1] * box.cpu_h[1] + box.cpu_h[4] * box.cpu_h[4] + box.cpu_h[7] * box.cpu_h[7]); + thickness[2] = + sqrt(box.cpu_h[2] * box.cpu_h[2] + box.cpu_h[5] * box.cpu_h[5] + box.cpu_h[8] * box.cpu_h[8]); int numCells[4]; numCells[0] = std::max(1, static_cast(ceil(thickness[0] * 0.2))); @@ -137,7 +141,8 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) double y12 = pos[n2 + N] - r1[1]; double z12 = pos[n2 + 2 * N] - r1[2]; applyMic(box, x12, y12, z12); - if (fabs(x12) > 2.0 || fabs(y12) > 2.0 || fabs(z12) > 2.0) continue; + if (fabs(x12) > 2.0 || fabs(y12) > 2.0 || fabs(z12) > 2.0) + continue; const double d2 = x12 * x12 + y12 * y12 + z12 * z12; if (d2 >= 4.0) continue; @@ -155,7 +160,7 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) } } - if (min_distance < 2.0) { + if (min_distance < 1.0) { printf( "Error: Minimum distance (%f Å) between atoms %d (%s) and %d (%s) is less than 1 Å.\n", min_distance, From 6337967a6ad6f2a28b70d1b922d3c9338d210717 Mon Sep 17 00:00:00 2001 From: tang070205 Date: Sat, 1 Feb 2025 12:41:43 +0800 Subject: [PATCH 10/20] Update check_distance.cu --- src/model/check_distance.cu | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index 98e63e81f..f7686c823 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -68,6 +68,14 @@ void findCell( cell[3] = cell[0] + numCells[0] * (cell[1] + numCells[1] * cell[2]); } +float getArea(const double* a, const double* b) +{ + const double s1 = a[1] * b[2] - a[2] * b[1]; + const double s2 = a[2] * b[0] - a[0] * b[2]; + const double s3 = a[0] * b[1] - a[1] * b[0]; + return sqrt(s1 * s1 + s2 * s2 + s3 * s3); +} + void calculate_min_atomic_distance(const Atom& atom, const Box& box) { const int N = atom.number_of_atoms; @@ -77,12 +85,16 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) int min_n1 = -1, min_n2 = -1; double thickness[3]; - thickness[0] = - sqrt(box.cpu_h[0] * box.cpu_h[0] + box.cpu_h[3] * box.cpu_h[3] + box.cpu_h[6] * box.cpu_h[6]); - thickness[1] = - sqrt(box.cpu_h[1] * box.cpu_h[1] + box.cpu_h[4] * box.cpu_h[4] + box.cpu_h[7] * box.cpu_h[7]); - thickness[2] = - sqrt(box.cpu_h[2] * box.cpu_h[2] + box.cpu_h[5] * box.cpu_h[5] + box.cpu_h[8] * box.cpu_h[8]); + double volume = abs( + box.cpu_h[0] * (box.cpu_h[4] * box.cpu_h[8] - box.cpu_h[5] * box.cpu_h[7]) + + box.cpu_h[1] * (box.cpu_h[5] * box.cpu_h[6] - box.cpu_h[3] * box.cpu_h[8]) + + box.cpu_h[2] * (box.cpu_h[3] * box.cpu_h[7] - box.cpu_h[4] * box.cpu_h[6])); + const double a[3] = {box.cpu_h[0], box.cpu_h[3], box.cpu_h[6]}; + const double b[3] = {box.cpu_h[1], box.cpu_h[4], box.cpu_h[7]}; + const double c[3] = {box.cpu_h[2], box.cpu_h[5], box.cpu_h[8]}; + thickness[0] = volume / getArea(b, c); + thickness[1] = volume / getArea(c, a); + thickness[2] = volume / getArea(a, b); int numCells[4]; numCells[0] = std::max(1, static_cast(ceil(thickness[0] * 0.2))); From ff8a9097f9bb10b42f37a5140e8a8115862b33a6 Mon Sep 17 00:00:00 2001 From: tang070205 Date: Sun, 2 Feb 2025 09:37:22 +0800 Subject: [PATCH 11/20] Update check_distance.cu --- src/model/check_distance.cu | 80 +++++++++++++------------------------ 1 file changed, 28 insertions(+), 52 deletions(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index f7686c823..d473993de 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -34,31 +34,25 @@ void applyMicOne(double& x12) void applyMic(const Box& box, double& x12, double& y12, double& z12) { - - double sx = box.cpu_h[9] * x12 + box.cpu_h[10] * y12 + box.cpu_h[11] * z12; - double sy = box.cpu_h[12] * x12 + box.cpu_h[13] * y12 + box.cpu_h[14] * z12; - double sz = box.cpu_h[15] * x12 + box.cpu_h[16] * y12 + box.cpu_h[17] * z12; - - if (box.pbc_x) - applyMicOne(sx); - if (box.pbc_y) - applyMicOne(sy); - if (box.pbc_z) - applyMicOne(sz); - - x12 = box.cpu_h[0] * sx + box.cpu_h[3] * sy + box.cpu_h[6] * sz; - y12 = box.cpu_h[1] * sx + box.cpu_h[4] * sy + box.cpu_h[7] * sz; - z12 = box.cpu_h[2] * sx + box.cpu_h[5] * sy + box.cpu_h[8] * sz; + int pbc[3] = {box.pbc_x, box.pbc_y, box.pbc_z}; + double s[3]; + for (int i = 0; i < 3; ++i) { + s[i] = box.cpu_h[9 + i * 3] * x12 + box.cpu_h[10 + i * 3] * y12 + box.cpu_h[11 + i * 3] * z12; + if (pbc[i]) + applyMicOne(s[i]); + } + x12 = box.cpu_h[0] * s[0] + box.cpu_h[3] * s[1] + box.cpu_h[6] * s[2]; + y12 = box.cpu_h[1] * s[0] + box.cpu_h[4] * s[1] + box.cpu_h[7] * s[2]; + z12 = box.cpu_h[2] * s[0] + box.cpu_h[5] * s[1] + box.cpu_h[8] * s[2]; } void findCell( const Box& box, const double* thickness, const double* r, const int* numCells, int* cell) { double s[3]; - s[0] = box.cpu_h[9] * r[0] + box.cpu_h[10] * r[1] + box.cpu_h[11] * r[2]; - s[1] = box.cpu_h[12] * r[0] + box.cpu_h[13] * r[1] + box.cpu_h[14] * r[2]; - s[2] = box.cpu_h[15] * r[0] + box.cpu_h[16] * r[1] + box.cpu_h[17] * r[2]; for (int d = 0; d < 3; ++d) { + s[d] = + box.cpu_h[9 + d * 3] * r[0] + box.cpu_h[10 + d * 3] * r[1] + box.cpu_h[11 + d * 3] * r[2]; cell[d] = floor(s[d] * thickness[d] * 0.2); if (cell[d] < 0) cell[d] += numCells[d]; @@ -68,14 +62,6 @@ void findCell( cell[3] = cell[0] + numCells[0] * (cell[1] + numCells[1] * cell[2]); } -float getArea(const double* a, const double* b) -{ - const double s1 = a[1] * b[2] - a[2] * b[1]; - const double s2 = a[2] * b[0] - a[0] * b[2]; - const double s3 = a[0] * b[1] - a[1] * b[0]; - return sqrt(s1 * s1 + s2 * s2 + s3 * s3); -} - void calculate_min_atomic_distance(const Atom& atom, const Box& box) { const int N = atom.number_of_atoms; @@ -84,28 +70,20 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) double min_distance = 5.0; int min_n1 = -1, min_n2 = -1; + int cell[4], numCells[4]; double thickness[3]; - double volume = abs( - box.cpu_h[0] * (box.cpu_h[4] * box.cpu_h[8] - box.cpu_h[5] * box.cpu_h[7]) + - box.cpu_h[1] * (box.cpu_h[5] * box.cpu_h[6] - box.cpu_h[3] * box.cpu_h[8]) + - box.cpu_h[2] * (box.cpu_h[3] * box.cpu_h[7] - box.cpu_h[4] * box.cpu_h[6])); - const double a[3] = {box.cpu_h[0], box.cpu_h[3], box.cpu_h[6]}; - const double b[3] = {box.cpu_h[1], box.cpu_h[4], box.cpu_h[7]}; - const double c[3] = {box.cpu_h[2], box.cpu_h[5], box.cpu_h[8]}; - thickness[0] = volume / getArea(b, c); - thickness[1] = volume / getArea(c, a); - thickness[2] = volume / getArea(a, b); - - int numCells[4]; - numCells[0] = std::max(1, static_cast(ceil(thickness[0] * 0.2))); - numCells[1] = std::max(1, static_cast(ceil(thickness[1] * 0.2))); - numCells[2] = std::max(1, static_cast(ceil(thickness[2] * 0.2))); + for (int i = 0; i < 3; ++i) { + thickness[i] = sqrt( + box.cpu_h[i] * box.cpu_h[i] + box.cpu_h[i + 3] * box.cpu_h[i + 3] + + box.cpu_h[i + 6] * box.cpu_h[i + 6]); + numCells[i] = std::max(1, static_cast(ceil(thickness[i] * 0.2))); + } numCells[3] = numCells[0] * numCells[1] * numCells[2]; - int cell[4]; - + std::vector cellContents(N, 0); std::vector cellCount(numCells[3], 0); std::vector cellCountSum(numCells[3], 0); + std::fill(cellCount.begin(), cellCount.end(), 0); for (int n = 0; n < N; ++n) { const double r[3] = {pos[n], pos[n + N], pos[n + 2 * N]}; @@ -117,9 +95,6 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) cellCountSum[i] = cellCountSum[i - 1] + cellCount[i - 1]; } - std::fill(cellCount.begin(), cellCount.end(), 0); - std::vector cellContents(N, 0); - for (int n = 0; n < N; ++n) { const double r[3] = {pos[n], pos[n + N], pos[n + 2 * N]}; findCell(box, thickness, r, numCells, cell); @@ -159,9 +134,9 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) if (d2 >= 4.0) continue; - double dist = sqrt(d2); - if (dist < min_distance) { - min_distance = dist; + double distance = d2; + if (distance < min_distance) { + min_distance = distance; min_n1 = n1; min_n2 = n2; } @@ -171,11 +146,12 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) } } } + double mini_distance = sqrt(min_distance); - if (min_distance < 1.0) { + if (mini_distance < 1.0) { printf( "Error: Minimum distance (%f Å) between atoms %d (%s) and %d (%s) is less than 1 Å.\n", - min_distance, + mini_distance, min_n1, atom.cpu_atom_symbol[min_n1].c_str(), min_n2, @@ -188,6 +164,6 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) atom.cpu_atom_symbol[min_n1].c_str(), min_n2, atom.cpu_atom_symbol[min_n2].c_str(), - min_distance); + mini_distance); } } From c67e19bc66f9a29aa6d37bb36c455341d62252c9 Mon Sep 17 00:00:00 2001 From: tang070205 Date: Sun, 2 Feb 2025 22:26:36 +0800 Subject: [PATCH 12/20] Update check_distance.cu --- src/model/check_distance.cu | 66 +++++++++++++++++++++---------------- 1 file changed, 37 insertions(+), 29 deletions(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index d473993de..29b9204e7 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -34,25 +34,31 @@ void applyMicOne(double& x12) void applyMic(const Box& box, double& x12, double& y12, double& z12) { - int pbc[3] = {box.pbc_x, box.pbc_y, box.pbc_z}; - double s[3]; - for (int i = 0; i < 3; ++i) { - s[i] = box.cpu_h[9 + i * 3] * x12 + box.cpu_h[10 + i * 3] * y12 + box.cpu_h[11 + i * 3] * z12; - if (pbc[i]) - applyMicOne(s[i]); - } - x12 = box.cpu_h[0] * s[0] + box.cpu_h[3] * s[1] + box.cpu_h[6] * s[2]; - y12 = box.cpu_h[1] * s[0] + box.cpu_h[4] * s[1] + box.cpu_h[7] * s[2]; - z12 = box.cpu_h[2] * s[0] + box.cpu_h[5] * s[1] + box.cpu_h[8] * s[2]; + + double sx = box.cpu_h[9] * x12 + box.cpu_h[10] * y12 + box.cpu_h[11] * z12; + double sy = box.cpu_h[12] * x12 + box.cpu_h[13] * y12 + box.cpu_h[14] * z12; + double sz = box.cpu_h[15] * x12 + box.cpu_h[16] * y12 + box.cpu_h[17] * z12; + + if (box.pbc_x) + applyMicOne(sx); + if (box.pbc_y) + applyMicOne(sy); + if (box.pbc_z) + applyMicOne(sz); + + x12 = box.cpu_h[0] * sx + box.cpu_h[3] * sy + box.cpu_h[6] * sz; + y12 = box.cpu_h[1] * sx + box.cpu_h[4] * sy + box.cpu_h[7] * sz; + z12 = box.cpu_h[2] * sx + box.cpu_h[5] * sy + box.cpu_h[8] * sz; } void findCell( const Box& box, const double* thickness, const double* r, const int* numCells, int* cell) { double s[3]; + s[0] = box.cpu_h[9] * r[0] + box.cpu_h[10] * r[1] + box.cpu_h[11] * r[2]; + s[1] = box.cpu_h[12] * r[0] + box.cpu_h[13] * r[1] + box.cpu_h[14] * r[2]; + s[2] = box.cpu_h[15] * r[0] + box.cpu_h[16] * r[1] + box.cpu_h[17] * r[2]; for (int d = 0; d < 3; ++d) { - s[d] = - box.cpu_h[9 + d * 3] * r[0] + box.cpu_h[10 + d * 3] * r[1] + box.cpu_h[11 + d * 3] * r[2]; cell[d] = floor(s[d] * thickness[d] * 0.2); if (cell[d] < 0) cell[d] += numCells[d]; @@ -67,17 +73,21 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) const int N = atom.number_of_atoms; const double* pos = atom.cpu_position_per_atom.data(); - double min_distance = 5.0; + double dist_sq = 5.0; int min_n1 = -1, min_n2 = -1; - int cell[4], numCells[4]; double thickness[3]; - for (int i = 0; i < 3; ++i) { - thickness[i] = sqrt( - box.cpu_h[i] * box.cpu_h[i] + box.cpu_h[i + 3] * box.cpu_h[i + 3] + - box.cpu_h[i + 6] * box.cpu_h[i + 6]); - numCells[i] = std::max(1, static_cast(ceil(thickness[i] * 0.2))); - } + thickness[0] = + sqrt(box.cpu_h[0] * box.cpu_h[0] + box.cpu_h[3] * box.cpu_h[3] + box.cpu_h[6] * box.cpu_h[6]); + thickness[1] = + sqrt(box.cpu_h[1] * box.cpu_h[1] + box.cpu_h[4] * box.cpu_h[4] + box.cpu_h[7] * box.cpu_h[7]); + thickness[2] = + sqrt(box.cpu_h[2] * box.cpu_h[2] + box.cpu_h[5] * box.cpu_h[5] + box.cpu_h[8] * box.cpu_h[8]); + + int cell[4], numCells[4]; + numCells[0] = std::max(1, static_cast(ceil(thickness[0] * 0.2))); + numCells[1] = std::max(1, static_cast(ceil(thickness[1] * 0.2))); + numCells[2] = std::max(1, static_cast(ceil(thickness[2] * 0.2))); numCells[3] = numCells[0] * numCells[1] * numCells[2]; std::vector cellContents(N, 0); @@ -134,11 +144,9 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) if (d2 >= 4.0) continue; - double distance = d2; - if (distance < min_distance) { - min_distance = distance; - min_n1 = n1; - min_n2 = n2; + double dist = d2; + if (dist < dist_sq) { + dist_sq = dist, min_n1 = n1, min_n2 = n2; } } } @@ -146,12 +154,12 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) } } } - double mini_distance = sqrt(min_distance); + double min_distance = sqrt(dist_sq); - if (mini_distance < 1.0) { + if (min_distance < 1.0) { printf( "Error: Minimum distance (%f Å) between atoms %d (%s) and %d (%s) is less than 1 Å.\n", - mini_distance, + min_distance, min_n1, atom.cpu_atom_symbol[min_n1].c_str(), min_n2, @@ -164,6 +172,6 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) atom.cpu_atom_symbol[min_n1].c_str(), min_n2, atom.cpu_atom_symbol[min_n2].c_str(), - mini_distance); + min_distance); } } From 6986d156a97a27b15c9dde60475ec0759d3d2032 Mon Sep 17 00:00:00 2001 From: tang070205 Date: Sun, 2 Feb 2025 22:39:49 +0800 Subject: [PATCH 13/20] Update check_distance.cu --- src/model/check_distance.cu | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index 29b9204e7..530d5c645 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -93,7 +93,6 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) std::vector cellContents(N, 0); std::vector cellCount(numCells[3], 0); std::vector cellCountSum(numCells[3], 0); - std::fill(cellCount.begin(), cellCount.end(), 0); for (int n = 0; n < N; ++n) { const double r[3] = {pos[n], pos[n + N], pos[n + 2 * N]}; @@ -104,7 +103,7 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) for (int i = 1; i < numCells[3]; ++i) { cellCountSum[i] = cellCountSum[i - 1] + cellCount[i - 1]; } - + std::fill(cellCount.begin(), cellCount.end(), 0); for (int n = 0; n < N; ++n) { const double r[3] = {pos[n], pos[n + N], pos[n + 2 * N]}; findCell(box, thickness, r, numCells, cell); From 5b2a371b5ab013d4b555bc85f67af443e4669612 Mon Sep 17 00:00:00 2001 From: tang070205 Date: Sun, 2 Feb 2025 22:41:32 +0800 Subject: [PATCH 14/20] Update check_distance.cu --- src/model/check_distance.cu | 66 ++++++++++++++++--------------------- 1 file changed, 29 insertions(+), 37 deletions(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index 530d5c645..86c584b90 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -34,31 +34,25 @@ void applyMicOne(double& x12) void applyMic(const Box& box, double& x12, double& y12, double& z12) { - - double sx = box.cpu_h[9] * x12 + box.cpu_h[10] * y12 + box.cpu_h[11] * z12; - double sy = box.cpu_h[12] * x12 + box.cpu_h[13] * y12 + box.cpu_h[14] * z12; - double sz = box.cpu_h[15] * x12 + box.cpu_h[16] * y12 + box.cpu_h[17] * z12; - - if (box.pbc_x) - applyMicOne(sx); - if (box.pbc_y) - applyMicOne(sy); - if (box.pbc_z) - applyMicOne(sz); - - x12 = box.cpu_h[0] * sx + box.cpu_h[3] * sy + box.cpu_h[6] * sz; - y12 = box.cpu_h[1] * sx + box.cpu_h[4] * sy + box.cpu_h[7] * sz; - z12 = box.cpu_h[2] * sx + box.cpu_h[5] * sy + box.cpu_h[8] * sz; + int pbc[3] = {box.pbc_x, box.pbc_y, box.pbc_z}; + double s[3]; + for (int i = 0; i < 3; ++i) { + s[i] = box.cpu_h[9 + i * 3] * x12 + box.cpu_h[10 + i * 3] * y12 + box.cpu_h[11 + i * 3] * z12; + if (pbc[i]) + applyMicOne(s[i]); + } + x12 = box.cpu_h[0] * s[0] + box.cpu_h[3] * s[1] + box.cpu_h[6] * s[2]; + y12 = box.cpu_h[1] * s[0] + box.cpu_h[4] * s[1] + box.cpu_h[7] * s[2]; + z12 = box.cpu_h[2] * s[0] + box.cpu_h[5] * s[1] + box.cpu_h[8] * s[2]; } void findCell( const Box& box, const double* thickness, const double* r, const int* numCells, int* cell) { double s[3]; - s[0] = box.cpu_h[9] * r[0] + box.cpu_h[10] * r[1] + box.cpu_h[11] * r[2]; - s[1] = box.cpu_h[12] * r[0] + box.cpu_h[13] * r[1] + box.cpu_h[14] * r[2]; - s[2] = box.cpu_h[15] * r[0] + box.cpu_h[16] * r[1] + box.cpu_h[17] * r[2]; for (int d = 0; d < 3; ++d) { + s[d] = + box.cpu_h[9 + d * 3] * r[0] + box.cpu_h[10 + d * 3] * r[1] + box.cpu_h[11 + d * 3] * r[2]; cell[d] = floor(s[d] * thickness[d] * 0.2); if (cell[d] < 0) cell[d] += numCells[d]; @@ -73,21 +67,17 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) const int N = atom.number_of_atoms; const double* pos = atom.cpu_position_per_atom.data(); - double dist_sq = 5.0; + double min_distance = 5.0; int min_n1 = -1, min_n2 = -1; - double thickness[3]; - thickness[0] = - sqrt(box.cpu_h[0] * box.cpu_h[0] + box.cpu_h[3] * box.cpu_h[3] + box.cpu_h[6] * box.cpu_h[6]); - thickness[1] = - sqrt(box.cpu_h[1] * box.cpu_h[1] + box.cpu_h[4] * box.cpu_h[4] + box.cpu_h[7] * box.cpu_h[7]); - thickness[2] = - sqrt(box.cpu_h[2] * box.cpu_h[2] + box.cpu_h[5] * box.cpu_h[5] + box.cpu_h[8] * box.cpu_h[8]); - int cell[4], numCells[4]; - numCells[0] = std::max(1, static_cast(ceil(thickness[0] * 0.2))); - numCells[1] = std::max(1, static_cast(ceil(thickness[1] * 0.2))); - numCells[2] = std::max(1, static_cast(ceil(thickness[2] * 0.2))); + double thickness[3]; + for (int i = 0; i < 3; ++i) { + thickness[i] = sqrt( + box.cpu_h[i] * box.cpu_h[i] + box.cpu_h[i + 3] * box.cpu_h[i + 3] + + box.cpu_h[i + 6] * box.cpu_h[i + 6]); + numCells[i] = std::max(1, static_cast(ceil(thickness[i] * 0.2))); + } numCells[3] = numCells[0] * numCells[1] * numCells[2]; std::vector cellContents(N, 0); @@ -143,9 +133,11 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) if (d2 >= 4.0) continue; - double dist = d2; - if (dist < dist_sq) { - dist_sq = dist, min_n1 = n1, min_n2 = n2; + double distance = d2; + if (distance < min_distance) { + min_distance = distance; + min_n1 = n1; + min_n2 = n2; } } } @@ -153,12 +145,12 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) } } } - double min_distance = sqrt(dist_sq); + double mini_distance = sqrt(min_distance); - if (min_distance < 1.0) { + if (mini_distance < 1.0) { printf( "Error: Minimum distance (%f Å) between atoms %d (%s) and %d (%s) is less than 1 Å.\n", - min_distance, + mini_distance, min_n1, atom.cpu_atom_symbol[min_n1].c_str(), min_n2, @@ -171,6 +163,6 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) atom.cpu_atom_symbol[min_n1].c_str(), min_n2, atom.cpu_atom_symbol[min_n2].c_str(), - min_distance); + mini_distance); } } From 7e0d2b4663813b288548147badade0cd2f18b652 Mon Sep 17 00:00:00 2001 From: tang070205 Date: Mon, 3 Feb 2025 18:24:52 +0800 Subject: [PATCH 15/20] Update check_distance.cu --- src/model/check_distance.cu | 57 +++++++++++++++++++++---------------- 1 file changed, 33 insertions(+), 24 deletions(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index 86c584b90..c9310b8b6 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -34,13 +34,17 @@ void applyMicOne(double& x12) void applyMic(const Box& box, double& x12, double& y12, double& z12) { - int pbc[3] = {box.pbc_x, box.pbc_y, box.pbc_z}; - double s[3]; - for (int i = 0; i < 3; ++i) { - s[i] = box.cpu_h[9 + i * 3] * x12 + box.cpu_h[10 + i * 3] * y12 + box.cpu_h[11 + i * 3] * z12; - if (pbc[i]) - applyMicOne(s[i]); - } + double sx = box.cpu_h[9] * x12 + box.cpu_h[10] * y12 + box.cpu_h[11] * z12; + double sy = box.cpu_h[12] * x12 + box.cpu_h[13] * y12 + box.cpu_h[14] * z12; + double sz = box.cpu_h[15] * x12 + box.cpu_h[16] * y12 + box.cpu_h[17] * z12; + + if (box.pbc_x) + applyMicOne(sx); + if (box.pbc_y) + applyMicOne(sy); + if (box.pbc_z) + applyMicOne(sz); + x12 = box.cpu_h[0] * s[0] + box.cpu_h[3] * s[1] + box.cpu_h[6] * s[2]; y12 = box.cpu_h[1] * s[0] + box.cpu_h[4] * s[1] + box.cpu_h[7] * s[2]; z12 = box.cpu_h[2] * s[0] + box.cpu_h[5] * s[1] + box.cpu_h[8] * s[2]; @@ -50,9 +54,10 @@ void findCell( const Box& box, const double* thickness, const double* r, const int* numCells, int* cell) { double s[3]; + s[0] = box.cpu_h[9] * r[0] + box.cpu_h[10] * r[1] + box.cpu_h[11] * r[2]; + s[1] = box.cpu_h[12] * r[0] + box.cpu_h[13] * r[1] + box.cpu_h[14] * r[2]; + s[2] = box.cpu_h[15] * r[0] + box.cpu_h[16] * r[1] + box.cpu_h[17] * r[2]; for (int d = 0; d < 3; ++d) { - s[d] = - box.cpu_h[9 + d * 3] * r[0] + box.cpu_h[10 + d * 3] * r[1] + box.cpu_h[11 + d * 3] * r[2]; cell[d] = floor(s[d] * thickness[d] * 0.2); if (cell[d] < 0) cell[d] += numCells[d]; @@ -67,17 +72,21 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) const int N = atom.number_of_atoms; const double* pos = atom.cpu_position_per_atom.data(); - double min_distance = 5.0; + double dist_sq = 5.0; int min_n1 = -1, min_n2 = -1; - int cell[4], numCells[4]; double thickness[3]; - for (int i = 0; i < 3; ++i) { - thickness[i] = sqrt( - box.cpu_h[i] * box.cpu_h[i] + box.cpu_h[i + 3] * box.cpu_h[i + 3] + - box.cpu_h[i + 6] * box.cpu_h[i + 6]); - numCells[i] = std::max(1, static_cast(ceil(thickness[i] * 0.2))); - } + thickness[0] = + sqrt(box.cpu_h[0] * box.cpu_h[0] + box.cpu_h[3] * box.cpu_h[3] + box.cpu_h[6] * box.cpu_h[6]); + thickness[1] = + sqrt(box.cpu_h[1] * box.cpu_h[1] + box.cpu_h[4] * box.cpu_h[4] + box.cpu_h[7] * box.cpu_h[7]); + thickness[2] = + sqrt(box.cpu_h[2] * box.cpu_h[2] + box.cpu_h[5] * box.cpu_h[5] + box.cpu_h[8] * box.cpu_h[8]); + + int cell[4], numCells[4]; + numCells[0] = std::max(1, static_cast(ceil(thickness[0] * 0.2))); + numCells[1] = std::max(1, static_cast(ceil(thickness[1] * 0.2))); + numCells[2] = std::max(1, static_cast(ceil(thickness[2] * 0.2))); numCells[3] = numCells[0] * numCells[1] * numCells[2]; std::vector cellContents(N, 0); @@ -133,9 +142,9 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) if (d2 >= 4.0) continue; - double distance = d2; - if (distance < min_distance) { - min_distance = distance; + double dist = d2; + if (dist < dist_sq) { + dist = dist_sq; min_n1 = n1; min_n2 = n2; } @@ -145,12 +154,12 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) } } } - double mini_distance = sqrt(min_distance); + double min_distance = sqrt(dist); - if (mini_distance < 1.0) { + if (min_distance < 1.0) { printf( "Error: Minimum distance (%f Å) between atoms %d (%s) and %d (%s) is less than 1 Å.\n", - mini_distance, + min_distance, min_n1, atom.cpu_atom_symbol[min_n1].c_str(), min_n2, @@ -163,6 +172,6 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) atom.cpu_atom_symbol[min_n1].c_str(), min_n2, atom.cpu_atom_symbol[min_n2].c_str(), - mini_distance); + min_distance); } } From e5551689de47bc2efb3b80186fd79e7bd8745cad Mon Sep 17 00:00:00 2001 From: tang070205 Date: Mon, 3 Feb 2025 18:31:56 +0800 Subject: [PATCH 16/20] Update check_distance.cu --- src/model/check_distance.cu | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index c9310b8b6..879bf82dd 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -45,9 +45,9 @@ void applyMic(const Box& box, double& x12, double& y12, double& z12) if (box.pbc_z) applyMicOne(sz); - x12 = box.cpu_h[0] * s[0] + box.cpu_h[3] * s[1] + box.cpu_h[6] * s[2]; - y12 = box.cpu_h[1] * s[0] + box.cpu_h[4] * s[1] + box.cpu_h[7] * s[2]; - z12 = box.cpu_h[2] * s[0] + box.cpu_h[5] * s[1] + box.cpu_h[8] * s[2]; + x12 = box.cpu_h[0] * sx + box.cpu_h[3] * sy + box.cpu_h[6] * sz; + y12 = box.cpu_h[1] * sx + box.cpu_h[4] * sy + box.cpu_h[7] * sz; + z12 = box.cpu_h[2] * sx + box.cpu_h[5] * sy + box.cpu_h[8] * sz; } void findCell( @@ -143,8 +143,8 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) continue; double dist = d2; - if (dist < dist_sq) { - dist = dist_sq; + if (dist_sq < dist) { + dist_sq = dist; min_n1 = n1; min_n2 = n2; } From b463a556a1b6ccedb5623cdb93f7e0fed67ebf5a Mon Sep 17 00:00:00 2001 From: tang070205 Date: Mon, 3 Feb 2025 18:32:40 +0800 Subject: [PATCH 17/20] Update check_distance.cu --- src/model/check_distance.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index 879bf82dd..f35db9dc0 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -154,7 +154,7 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) } } } - double min_distance = sqrt(dist); + double min_distance = sqrt(dist_sq); if (min_distance < 1.0) { printf( From bf27c5e316e8fa009f87ee9e81c31b3bf1fc4f3b Mon Sep 17 00:00:00 2001 From: tang070205 Date: Mon, 3 Feb 2025 18:48:52 +0800 Subject: [PATCH 18/20] Update check_distance.cu --- src/model/check_distance.cu | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index f35db9dc0..261a35b6c 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -34,6 +34,7 @@ void applyMicOne(double& x12) void applyMic(const Box& box, double& x12, double& y12, double& z12) { + double sx = box.cpu_h[9] * x12 + box.cpu_h[10] * y12 + box.cpu_h[11] * z12; double sy = box.cpu_h[12] * x12 + box.cpu_h[13] * y12 + box.cpu_h[14] * z12; double sz = box.cpu_h[15] * x12 + box.cpu_h[16] * y12 + box.cpu_h[17] * z12; @@ -143,7 +144,7 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) continue; double dist = d2; - if (dist_sq < dist) { + if (dist < dist_sq) { dist_sq = dist; min_n1 = n1; min_n2 = n2; From a5309c519e59824eabc264523afd20adb2249180 Mon Sep 17 00:00:00 2001 From: tang070205 Date: Mon, 3 Feb 2025 21:02:37 +0800 Subject: [PATCH 19/20] modify the real coordinates of the applyMic function --- src/model/check_distance.cu | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index 261a35b6c..75fddc4ef 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -46,9 +46,9 @@ void applyMic(const Box& box, double& x12, double& y12, double& z12) if (box.pbc_z) applyMicOne(sz); - x12 = box.cpu_h[0] * sx + box.cpu_h[3] * sy + box.cpu_h[6] * sz; - y12 = box.cpu_h[1] * sx + box.cpu_h[4] * sy + box.cpu_h[7] * sz; - z12 = box.cpu_h[2] * sx + box.cpu_h[5] * sy + box.cpu_h[8] * sz; + x12 = box.cpu_h[0] * sx + box.cpu_h[1] * sy + box.cpu_h[2] * sz; + y12 = box.cpu_h[3] * sx + box.cpu_h[4] * sy + box.cpu_h[5] * sz; + z12 = box.cpu_h[6] * sx + box.cpu_h[7] * sy + box.cpu_h[8] * sz; } void findCell( @@ -166,7 +166,7 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) min_n2, atom.cpu_atom_symbol[min_n2].c_str()); PRINT_INPUT_ERROR("There are two atoms with a distance less than 1 Å."); - } else if (min_n1 != -1 && min_n2 != -1) { + } else { printf( "Minimum distance between atoms %d (%s) and %d (%s): %f Å\n", min_n1, From aca7d333c0ff4ead01073d4a2f70f53dceab8475 Mon Sep 17 00:00:00 2001 From: tang070205 Date: Mon, 3 Feb 2025 21:41:00 +0800 Subject: [PATCH 20/20] change to only check the minimum distance without report errors --- src/model/check_distance.cu | 26 ++++++++------------------ 1 file changed, 8 insertions(+), 18 deletions(-) diff --git a/src/model/check_distance.cu b/src/model/check_distance.cu index 75fddc4ef..c284daf5f 100644 --- a/src/model/check_distance.cu +++ b/src/model/check_distance.cu @@ -157,22 +157,12 @@ void calculate_min_atomic_distance(const Atom& atom, const Box& box) } double min_distance = sqrt(dist_sq); - if (min_distance < 1.0) { - printf( - "Error: Minimum distance (%f Å) between atoms %d (%s) and %d (%s) is less than 1 Å.\n", - min_distance, - min_n1, - atom.cpu_atom_symbol[min_n1].c_str(), - min_n2, - atom.cpu_atom_symbol[min_n2].c_str()); - PRINT_INPUT_ERROR("There are two atoms with a distance less than 1 Å."); - } else { - printf( - "Minimum distance between atoms %d (%s) and %d (%s): %f Å\n", - min_n1, - atom.cpu_atom_symbol[min_n1].c_str(), - min_n2, - atom.cpu_atom_symbol[min_n2].c_str(), - min_distance); - } + printf( + "Minimum distance between atoms %d (%s) and %d (%s): %f Å\n", + min_n1, + atom.cpu_atom_symbol[min_n1].c_str(), + min_n2, + atom.cpu_atom_symbol[min_n2].c_str(), + min_distance); + }