From 647380a3578ca49c95675c6c1437a838219351e0 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Wed, 20 Oct 2021 13:21:07 -0600 Subject: [PATCH] Avoid bitshift that gave incorrect results on GPU --- src/KOKKOS/compute_orientorder_atom_kokkos.cpp | 12 +++++++----- src/compute_orientorder_atom.cpp | 8 +++++--- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/KOKKOS/compute_orientorder_atom_kokkos.cpp b/src/KOKKOS/compute_orientorder_atom_kokkos.cpp index 6edb500870..548c503628 100644 --- a/src/KOKKOS/compute_orientorder_atom_kokkos.cpp +++ b/src/KOKKOS/compute_orientorder_atom_kokkos.cpp @@ -26,7 +26,7 @@ #include "atom_masks.h" #include "kokkos.h" #include "math_const.h" -#include "math_special.h" +#include "math_special_kokkos.h" #include "memory_kokkos.h" #include "neigh_list.h" #include "neigh_request.h" @@ -38,7 +38,7 @@ using namespace LAMMPS_NS; using namespace MathConst; -using namespace MathSpecial; +using namespace MathSpecialKokkos; using namespace std; #ifdef DBL_EPSILON @@ -556,8 +556,9 @@ void ComputeOrientOrderAtomKokkos::calc_boop2(int ncount, int ii) co // structure enforces visiting only one member of each // such symmetry (invariance) group. - const int sgn = 1 - 2*(m1&1); - // m1 <= 0, and Qlm[-m] = (-1)^m*conjg(Qlm[m]). sgn = (-1)^m. + //const int sgn = 1 - 2*(m1&1); + const int sgn = powint(-1,m1); // sgn = (-1)^m + // m1 <= 0, and Qlm[-m] = (-1)^m*conjg(Qlm[m]) SNAcomplex Q1Q2; Q1Q2.re = (d_qnm(ii,il,-m1).re*d_qnm(ii,il,m2).re + d_qnm(ii,il,-m1).im*d_qnm(ii,il,m2).im)*sgn; Q1Q2.im = (d_qnm(ii,il,-m1).re*d_qnm(ii,il,m2).im - d_qnm(ii,il,-m1).im*d_qnm(ii,il,m2).re)*sgn; @@ -601,7 +602,8 @@ void ComputeOrientOrderAtomKokkos::calc_boop2(int ncount, int ii) co for (int m = -l; m < 0; m++) { // Computed only qnm for m>=0. // qnm[-m] = (-1)^m * conjg(qnm[m]) - const int sgn = 1 - 2*(m&1); // sgn = (-1)^m + //const int sgn = 1 - 2*(m&1); // sgn = (-1)^m + const int sgn = powint(-1,m); // sgn = (-1)^m d_qnarray(i,jj++) = d_qnm(ii,il,-m).re * qnfac * sgn; d_qnarray(i,jj++) = -d_qnm(ii,il,-m).im * qnfac * sgn; } diff --git a/src/compute_orientorder_atom.cpp b/src/compute_orientorder_atom.cpp index 4cca0316c5..876ee8e3ed 100644 --- a/src/compute_orientorder_atom.cpp +++ b/src/compute_orientorder_atom.cpp @@ -532,8 +532,9 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, // structure enforces visiting only one member of each // such symmetry (invariance) group. - const int sgn = 1 - 2*(m1&1); - // m1 <= 0, and Qlm[-m] = (-1)^m*conjg(Qlm[m]). sgn = (-1)^m. + //const int sgn = 1 - 2*(m1&1); + const int sgn = powint(-1,m1); // sgn = (-1)^m + // m1 <= 0, and Qlm[-m] = (-1)^m*conjg(Qlm[m]) const double Q1Q2_r = (qnm_r[il][-m1]*qnm_r[il][m2] + qnm_i[il][-m1]*qnm_i[il][m2])*sgn; const double Q1Q2_i = (qnm_r[il][-m1]*qnm_i[il][m2] - qnm_i[il][-m1]*qnm_r[il][m2])*sgn; const double Q1Q2Q3 = Q1Q2_r*qnm_r[il][m3] - Q1Q2_i*qnm_i[il][m3]; @@ -576,7 +577,8 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, for (int m = -l; m < 0; m++) { // Computed only qnm for m>=0. // qnm[-m] = (-1)^m * conjg(qnm[m]) - const int sgn = 1 - 2*(m&1); // sgn = (-1)^m + //const int sgn = 1 - 2*(m&1); // sgn = (-1)^m + const int sgn = powint(-1,m); // sgn = (-1)^m qn[jj++] = qnm_r[il][-m] * qnfac * sgn; qn[jj++] = -qnm_i[il][-m] * qnfac * sgn; }