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

Use lighter construct than powint

Browse Source
This commit is contained in:
Stan Gerald Moore
2021-10-27 13:41:29 -06:00
parent 2d52821cf5
commit 2ca7dcb853
2 changed files with 14 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
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_kokkos.h" #include "math_special.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 MathSpecialKokkos; using namespace MathSpecial;
using namespace std; using namespace std;
#ifdef DBL_EPSILON #ifdef DBL_EPSILON
@ -546,6 +546,7 @@ void ComputeOrientOrderAtomKokkos<DeviceType>::calc_boop2(int ncount, int ii) co
for (int il = 0; il < nqlist; il++) { for (int il = 0; il < nqlist; il++) {
int l = d_qlist[il]; int l = d_qlist[il];
double wlsum = 0.0; double wlsum = 0.0;
int sgn = 1;
for (int m1 = -l; m1 <= 0; m1++) { for (int m1 = -l; m1 <= 0; m1++) {
for (int m2 = 0; m2 <= ((-m1)>>1); m2++) { for (int m2 = 0; m2 <= ((-m1)>>1); m2++) {
const int m3 = -(m1 + m2); const int m3 = -(m1 + m2);
@ -556,8 +557,6 @@ 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 = powint(-1,m1); // sgn = (-1)^m
// m1 <= 0, and Qlm[-m] = (-1)^m*conjg(Qlm[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;
@ -565,6 +564,8 @@ void ComputeOrientOrderAtomKokkos<DeviceType>::calc_boop2(int ncount, int ii) co
const double Q1Q2Q3 = Q1Q2.re*d_qnm(ii,il,m3).re - Q1Q2.im*d_qnm(ii,il,m3).im; const double Q1Q2Q3 = Q1Q2.re*d_qnm(ii,il,m3).re - Q1Q2.im*d_qnm(ii,il,m3).im;
const double c = d_w3jlist[widx_count++]; const double c = d_w3jlist[widx_count++];
wlsum += Q1Q2Q3*c; wlsum += Q1Q2Q3*c;
sgn = -sgn; // sgn = (-1)^m
} }
} }
d_qnarray(i,jj++) = wlsum/d_qnormfac2(il); d_qnarray(i,jj++) = wlsum/d_qnormfac2(il);
@ -599,13 +600,14 @@ void ComputeOrientOrderAtomKokkos<DeviceType>::calc_boop2(int ncount, int ii) co
} }
else { else {
const double qnfac = d_qnormfac(il)/d_qnarray(i,il); const double qnfac = d_qnormfac(il)/d_qnarray(i,il);
int sgn = 1;
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 = 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;
sgn = -sgn; // sgn = (-1)^m
} }
for (int m = 0; m < l+1; m++) { for (int m = 0; m < l+1; m++) {
d_qnarray(i,jj++) = d_qnm(ii,il,m).re * qnfac; d_qnarray(i,jj++) = d_qnm(ii,il,m).re * qnfac;

10
src/compute_orientorder_atom.cpp
View File

@ -524,6 +524,7 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
for (int il = 0; il < nqlist; il++) { for (int il = 0; il < nqlist; il++) {
int l = qlist[il]; int l = qlist[il];
double wlsum = 0.0; double wlsum = 0.0;
int sgn = 1;
for (int m1 = -l; m1 <= 0; m1++) { for (int m1 = -l; m1 <= 0; m1++) {
for (int m2 = 0; m2 <= ((-m1)>>1); m2++) { for (int m2 = 0; m2 <= ((-m1)>>1); m2++) {
const int m3 = -(m1 + m2); const int m3 = -(m1 + m2);
@ -534,14 +535,14 @@ 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 = powint(-1,m1); // sgn = (-1)^m
// m1 <= 0, and Qlm[-m] = (-1)^m*conjg(Qlm[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];
const double c = w3jlist[widx_count++]; const double c = w3jlist[widx_count++];
wlsum += Q1Q2Q3*c; wlsum += Q1Q2Q3*c;
sgn = -sgn; // sgn = (-1)^m
} }
} }
qn[jj++] = wlsum/qnormfac2[il]; qn[jj++] = wlsum/qnormfac2[il];
@ -576,13 +577,14 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
} }
else { else {
double qnfac = qnormfac[il]/qn[il]; double qnfac = qnormfac[il]/qn[il];
int sgn = 1;
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 = 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;
sgn = -sgn; // sgn = (-1)^m
} }
for (int m = 0; m < l+1; m++) { for (int m = 0; m < l+1; m++) {
qn[jj++] = qnm_r[il][m] * qnfac; qn[jj++] = qnm_r[il][m] * qnfac;