ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Avoid bitshift that gave incorrect results on GPU

Browse Source
This commit is contained in:
Stan Gerald Moore
2021-10-20 13:21:07 -06:00
parent 01d96fc684
commit 647380a357
2 changed files with 12 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
src/KOKKOS/compute_orientorder_atom_kokkos.cpp
View File

@ -26,7 +26,7 @@
#include "atom_masks.h" #include "atom_masks.h"
#include "kokkos.h" #include "kokkos.h"
#include "math_const.h" #include "math_const.h"
#include "math_special.h" #include "math_special_kokkos.h"
#include "memory_kokkos.h" #include "memory_kokkos.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "neigh_request.h" #include "neigh_request.h"
@ -38,7 +38,7 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst; using namespace MathConst;
using namespace MathSpecial; using namespace MathSpecialKokkos;
using namespace std; using namespace std;
#ifdef DBL_EPSILON #ifdef DBL_EPSILON
@ -556,8 +556,9 @@ void ComputeOrientOrderAtomKokkos<DeviceType>::calc_boop2(int ncount, int ii) co
// structure enforces visiting only one member of each // structure enforces visiting only one member of each
// such symmetry (invariance) group. // such symmetry (invariance) group.
const int sgn = 1 - 2*(m1&1); //const int sgn = 1 - 2*(m1&1);
// m1 <= 0, and Qlm[-m] = (-1)^m*conjg(Qlm[m]). sgn = (-1)^m. const int sgn = powint(-1,m1); // sgn = (-1)^m
// m1 <= 0, and Qlm[-m] = (-1)^m*conjg(Qlm[m])
SNAcomplex Q1Q2; 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.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; 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<DeviceType>::calc_boop2(int ncount, int ii) co
for (int m = -l; m < 0; m++) { for (int m = -l; m < 0; m++) {
// Computed only qnm for m>=0. // Computed only qnm for m>=0.
// qnm[-m] = (-1)^m * conjg(qnm[m]) // 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).re * qnfac * sgn;
d_qnarray(i,jj++) = -d_qnm(ii,il,-m).im * qnfac * sgn; d_qnarray(i,jj++) = -d_qnm(ii,il,-m).im * qnfac * sgn;
} }

8
src/compute_orientorder_atom.cpp
View File

@ -532,8 +532,9 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
// structure enforces visiting only one member of each // structure enforces visiting only one member of each
// such symmetry (invariance) group. // such symmetry (invariance) group.
const int sgn = 1 - 2*(m1&1); //const int sgn = 1 - 2*(m1&1);
// m1 <= 0, and Qlm[-m] = (-1)^m*conjg(Qlm[m]). sgn = (-1)^m. 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_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 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]; 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++) { for (int m = -l; m < 0; m++) {
// Computed only qnm for m>=0. // Computed only qnm for m>=0.
// qnm[-m] = (-1)^m * conjg(qnm[m]) // 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_r[il][-m] * qnfac * sgn;
qn[jj++] = -qnm_i[il][-m] * qnfac * sgn; qn[jj++] = -qnm_i[il][-m] * qnfac * sgn;
} }