55 int unique = (gbisParams->
numPatches == 1) ? 1 : 0;
59 int cid = gbisParams->
cid;
65 float a_cut = gbisParams->
a_cut - fsMax;
66 float a_cut_ps2 = (a_cut+fsMax)*(a_cut+fsMax);
67 float r_cut = gbisParams->
cutoff;
68 float a_cut2 = a_cut*a_cut;
69 float r_cut2 = r_cut*r_cut;
73 double t1 = 1.0*clock()/CLOCKS_PER_SEC;
79 float rhoi0, rhois, rhoj0, rhojs, rhois2, rhojs2, rhois2_16;
80 float ami2, amj2, api2, apj2;
81 int numGBISPairlists = 4;
83 float max_gcut2 = r_cut;
84 if (a_cut+fsMax > r_cut)
85 max_gcut2 = a_cut+fsMax;
88 max_gcut2 *= max_gcut2;
90 int maxGroupPairs = params->
numAtoms[1];
91 short *groupPairs =
new short[maxGroupPairs];
94 for (
int ngi = minIg; ngi < maxI; ) {
95 int numGroupPairs = 0;
105 dr = ngri.
x - ngrj.
x;
107 dr = ngri.
y - ngrj.
y;
109 dr = ngri.
z - ngrj.
z;
111 if (r2 < max_gcut2) {
112 groupPairs[numGroupPairs++] = ngj;
118 for (
int i=ngi; i<ngi+iGroupSize; i++) {
121 int *size =
new int[numGBISPairlists];
123 for (
int k = 0; k < numGBISPairlists; k++) {
134 rhois = gbisParams->
intRad[0][2*i+1];
135 api2 = a_cut + rhois;
137 ami2 = a_cut - rhois;
139 rhois2 = rhois*rhois;
140 rhois2_16 = 16.0*rhois2;
144 for (
int j = i+1; j < ngi+iGroupSize; j++) {
157 rhojs = gbisParams->
intRad[1][2*j+1];
158 rhojs2 = rhojs*rhojs;
160 amj2 = a_cut - rhojs;
163 apj2 = a_cut + rhojs;
167 if ( r2 < ami2 && r2 < amj2 &&
168 r2 > rhois2_16 && r2 > 16.0*rhojs2 ) {
170 pairs[0][size[0]++] = j;
171 }
else if ( r2 < api2 && r2 < apj2 &&
172 r2 > ami2 && r2 > amj2 ) {
174 pairs[1][size[1]++] = j;
175 }
else if ( r2 < a_cut_ps2 ) {
177 pairs[2][size[2]++] = j;
178 }
else if ( r2 < r_cut2 ) {
180 pairs[3][size[3]++] = j;
185 for (
int g = 0; g < numGroupPairs; g++) {
186 int ngj = groupPairs[g];
189 int maxJ = ngj+jGroupSize;
190 for (
int j = ngj; j < maxJ; j++) {
205 rhojs = gbisParams->
intRad[1][2*j+1];
206 rhojs2 = rhojs*rhojs;
208 amj2 = a_cut - rhojs;
211 apj2 = a_cut + rhojs;
214 if ( r2 < ami2 && r2 < amj2 &&
215 r2 > rhois2_16 && r2 > 16.0*rhojs2 ) {
218 pairs[0][size[0]++] = j;
219 }
else if ( r2 < api2 && r2 < apj2 &&
220 r2 > ami2 && r2 > amj2 ) {
223 pairs[1][size[1]++] = j;
224 }
else if ( r2 < a_cut_ps2 ) {
227 pairs[2][size[2]++] = j;
228 }
else if ( r2 < r_cut2 ) {
231 pairs[3][size[3]++] = j;
235 for (
int k = 0; k < numGBISPairlists; k++) {
244 for (
int s = 0; s < strideIg; s++) {
249 for (
int k = 0; k < numGBISPairlists; k++)
253 double t2 = 1.0*clock()/CLOCKS_PER_SEC;
254 CkPrintf(
"PHASELST: %8.3f ms @ %8.3f ns/iter for %i iter\n",1000.0*(t2-t1), 1000000000.0*(t2-t1)/nops, nops);
Pairlists * gbisStepPairlists[4]
void newsize(int list_size)