ComputeGBIS.C File Reference

#include "ComputeNonbondedUtil.h"
#include "ComputeGBIS.inl"
#include "time.h"

Go to the source code of this file.

Defines

#define CHECK_PRIORITIES   0

Functions

void pairlistFromAll (nonbonded *params, GBISParamStruct *gbisParams, int minIg, int strideIg, int maxI)


Define Documentation

#define CHECK_PRIORITIES   0

Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved.

Definition at line 9 of file ComputeGBIS.C.


Function Documentation

void pairlistFromAll ( nonbonded params,
GBISParamStruct gbisParams,
int  minIg,
int  strideIg,
int  maxI 
) [inline]

Definition at line 44 of file ComputeGBIS.C.

References GBISParamStruct::a_cut, GBISParamStruct::cid, GBISParamStruct::cutoff, FS_MAX, GBISParamStruct::gbisPhase, GBISParamStruct::gbisStepPairlists, if(), GBISParamStruct::intRad, j, GBISParamStruct::maxGroupRadius, Pairlists::newsize(), CompAtom::nonbondedGroupSize, nonbonded::numAtoms, GBISParamStruct::numPatches, nonbonded::offset, nonbonded::p, CompAtom::position, GBISParamStruct::sequence, Vector::x, Vector::y, and Vector::z.

Referenced by ComputeNonbondedUtil::calcGBIS().

00050   {
00051   const BigReal offset_x = params->offset.x;
00052   const BigReal offset_y = params->offset.y;
00053   const BigReal offset_z = params->offset.z;
00054 
00055   int unique = (gbisParams->numPatches == 1) ? 1 : 0;
00056 
00057   int maxPairs = params->numAtoms[1];
00058   int seq = gbisParams->sequence;
00059   int cid = gbisParams->cid;
00060   int pe = CkMyPe();
00061   int gbisPhase = gbisParams->gbisPhase;
00062 
00063   //determine long and short cutoff
00064   float fsMax = FS_MAX; // gbisParams->fsMax;
00065   float a_cut = gbisParams->a_cut - fsMax;
00066   float a_cut_ps2 = (a_cut+fsMax)*(a_cut+fsMax);
00067   float r_cut = gbisParams->cutoff;
00068   float a_cut2 = a_cut*a_cut;
00069   float r_cut2 = r_cut*r_cut;
00070   int s = 0;// 0 is short cutoff list
00071   int l = 1;// 1 is long-short cutoff list
00072 #ifdef BENCHMARK
00073   double t1 = 1.0*clock()/CLOCKS_PER_SEC;
00074   int nops = 0;
00075 #endif
00076   Position ri, rj;
00077   Position ngri, ngrj;
00078   float r, dr, r2;
00079   float rhoi0, rhois, rhoj0, rhojs, rhois2, rhojs2, rhois2_16;
00080   float ami2, amj2, api2, apj2;
00081   int numGBISPairlists = 4;
00082 
00083   float max_gcut2  = r_cut;
00084   if (a_cut+fsMax > r_cut)
00085     max_gcut2 = a_cut+fsMax;
00086   max_gcut2 += 2.0*gbisParams->maxGroupRadius;
00087   //CkPrintf("max_cut = %f, %f\n",max_gcut2,gbisParams->maxGroupRadius);
00088   max_gcut2 *= max_gcut2;
00089 
00090   int maxGroupPairs = params->numAtoms[1];
00091   short *groupPairs = new short[maxGroupPairs];/*delete*/
00092 
00093   //foreach nonbonded group i
00094   for (int ngi = minIg; ngi < maxI; /*ngi updated at loop bottom*/ ) {
00095     int numGroupPairs = 0;
00096     ngri = params->p[0][ngi].position;
00097     ngri.x += offset_x;
00098     ngri.y += offset_y;
00099     ngri.z += offset_z;
00100 
00101     //find close j-groups; include i-group
00102     for (int ngj = unique*(ngi+params->p[0][ngi].nonbondedGroupSize); ngj < params->numAtoms[1]; ngj+=params->p[1][ngj].nonbondedGroupSize) {
00103       
00104       ngrj = params->p[1][ngj].position;
00105       dr = ngri.x - ngrj.x;
00106       r2 = dr*dr;
00107       dr = ngri.y - ngrj.y;
00108       r2 += dr*dr;
00109       dr = ngri.z - ngrj.z;
00110       r2 += dr*dr;
00111       if (r2 < max_gcut2) {
00112         groupPairs[numGroupPairs++] = ngj;
00113       }
00114     }
00115 
00116     //for all i in i-group
00117     int iGroupSize = params->p[0][ngi].nonbondedGroupSize;
00118     for (int i=ngi; i<ngi+iGroupSize; i++) {
00119       //CkPrintf("\tFORALL i=%05i\n",params->pExt[0][i].id);
00120       //extend pairlists
00121       int *size = new int[numGBISPairlists];
00122       plint **pairs = new plint*[numGBISPairlists];
00123       for (int k = 0; k < numGBISPairlists; k++) {
00124         size[k] = 0;
00125         pairs[k] = gbisParams->gbisStepPairlists[k]->
00126                              newlist(maxPairs);
00127       }
00128 
00129       //load i values
00130       ri = params->p[0][i].position;
00131       ri.x += offset_x;
00132       ri.y += offset_y;
00133       ri.z += offset_z;
00134       rhois = gbisParams->intRad[0][2*i+1];
00135       api2 = a_cut + rhois;
00136       api2 *= api2;
00137       ami2 = a_cut - rhois;
00138       ami2 *= ami2;
00139       rhois2 = rhois*rhois;
00140       rhois2_16 = 16.0*rhois2;
00141 
00142       //for all other i's in i-group
00143       if (unique)
00144       for (int j = i+1; j < ngi+iGroupSize; j++) {
00145 
00146 #ifdef BENCHMARK
00147           nops ++;
00148 #endif
00149           rj = params->p[1][j].position;
00150           dr = ri.x - rj.x;
00151           r2 = dr*dr;
00152           dr = ri.y - rj.y;
00153           r2 += dr*dr;
00154           dr = ri.z - rj.z;
00155           r2 += dr*dr;
00156 
00157           rhojs = gbisParams->intRad[1][2*j+1];
00158           rhojs2 = rhojs*rhojs;
00159           //r = sqrt(r2);
00160           amj2 = a_cut - rhojs;
00161           amj2 *= amj2;
00162 
00163           apj2 = a_cut + rhojs;
00164           apj2 *= apj2;
00165           //CkPrintf("IPAIR %5i %5i %5i\n",0*gbisParams->cid,params->pExt[0][i].id,params->pExt[1][j].id);
00166 
00167           if (  r2 < ami2 && r2 < amj2 &&
00168             r2 > rhois2_16 && r2 > 16.0*rhojs2 ) {
00169             //belongs to 22
00170             pairs[0][size[0]++] = j;
00171           } else if ( r2 < api2 && r2 < apj2 &&
00172             r2 > ami2 && r2 > amj2 ) {
00173             //belongs to 11
00174             pairs[1][size[1]++] = j;
00175           } else if ( r2 < a_cut_ps2 ) {
00176             //belongs to other a_cut
00177             pairs[2][size[2]++] = j;
00178           } else if ( r2 < r_cut2 ) {
00179             //belongs to r_cut and (r_cut > a_cut)
00180             pairs[3][size[3]++] = j;
00181           }
00182       }
00183 
00184       //foreach j-group
00185       for (int g = 0; g < numGroupPairs; g++) {
00186         int ngj = groupPairs[g];
00187         int jGroupSize = params->p[1][ngj].nonbondedGroupSize;
00188         //for each j in j-group
00189         int maxJ = ngj+jGroupSize;
00190         for (int j = ngj; j < maxJ; j++) {
00191 
00192           //CkPrintf("for j-atom%i\n",j);
00193 #ifdef BENCHMARK
00194           nops ++;
00195 #endif
00196           rj = params->p[1][j].position;
00197           dr = ri.x - rj.x;
00198           r2 = dr*dr;
00199           dr = ri.y - rj.y;
00200           r2 += dr*dr;
00201           dr = ri.z - rj.z;
00202           r2 += dr*dr;
00203           //CkPrintf("PAIR %5i %5i %5i\n",0*gbisParams->cid,params->pExt[0][i].id,params->pExt[1][j].id);
00204 
00205           rhojs = gbisParams->intRad[1][2*j+1];
00206           rhojs2 = rhojs*rhojs;
00207           //r = sqrt(r2);
00208           amj2 = a_cut - rhojs;
00209           amj2 *= amj2;
00210 
00211           apj2 = a_cut + rhojs;
00212           apj2 *= apj2;
00213 
00214           if (  r2 < ami2 && r2 < amj2 &&
00215             r2 > rhois2_16 && r2 > 16.0*rhojs2 ) {
00216             //belongs to 22
00217             //CkPrintf("Pair22 %06i %06i\n",params->pExt[0][i].id,params->pExt[1][j].id);
00218             pairs[0][size[0]++] = j;
00219           } else if ( r2 < api2 && r2 < apj2 &&
00220             r2 > ami2 && r2 > amj2 ) {
00221             //CkPrintf("Pair11 %06i %06i\n",params->pExt[0][i].id,params->pExt[1][j].id);
00222             //belongs to 11
00223             pairs[1][size[1]++] = j;
00224           } else if ( r2 < a_cut_ps2 ) {
00225             //CkPrintf("PairAA %06i %06i\n",params->pExt[0][i].id,params->pExt[1][j].id);
00226             //belongs to other a_cut
00227             pairs[2][size[2]++] = j;
00228           } else if ( r2 < r_cut2 ) {
00229             //CkPrintf("PairRR %06i %06i\n",params->pExt[0][i].id,params->pExt[1][j].id);
00230             //belongs to r_cut and (r_cut > a_cut)
00231             pairs[3][size[3]++] = j;
00232           }
00233         }//end j inner loop
00234       }//end j-group loop
00235     for (int k = 0; k < numGBISPairlists; k++) {
00236       gbisParams->gbisStepPairlists[k]->newsize(size[k]);
00237     }
00238     delete[] size;
00239     delete[] pairs;
00240     }//end i all atom loop
00241 
00242     //jump to next nbg for round-robin
00243     //same iteration patter followed in all three gbisPhases below
00244     for (int s = 0; s < strideIg; s++) {
00245       ngi+=params->p[0][ngi].nonbondedGroupSize;
00246     }
00247   }//end i-group loop
00248     delete[] groupPairs;
00249   for (int k = 0; k < numGBISPairlists; k++)
00250     gbisParams->gbisStepPairlists[k]->reset();
00251 
00252 #ifdef BENCHMARK
00253   double t2 = 1.0*clock()/CLOCKS_PER_SEC;
00254   CkPrintf("PHASELST: %8.3f ms @ %8.3f ns/iter for %i iter\n",1000.0*(t2-t1), 1000000000.0*(t2-t1)/nops, nops);
00255 #endif
00256 }//end pairlist method


Generated on Sat Nov 18 01:17:16 2017 for NAMD by  doxygen 1.4.7