ComputeGBIS.C File Reference

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

Go to the source code of this file.

Macros
#define	CHECK_PRIORITIES   0

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

Macro Definition 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, GBISParamStruct::intRad, GBISParamStruct::maxGroupRadius, Pairlists::newsize(), CompAtom::nonbondedGroupSize, nonbonded::numAtoms, GBISParamStruct::numPatches, nonbonded::offset, nonbonded::p, CompAtom::position, Pairlists::reset(), GBISParamStruct::sequence, Vector::x, Vector::y, and Vector::z.

Referenced by ComputeNonbondedUtil::calcGBIS().

  {
  const BigReal offset_x = params->offset.x;
  const BigReal offset_y = params->offset.y;
  const BigReal offset_z = params->offset.z;

  int unique = (gbisParams->numPatches == 1) ? 1 : 0;

  int maxPairs = params->numAtoms[1];
  int seq = gbisParams->sequence;
  int cid = gbisParams->cid;
  int pe = CkMyPe();
  int gbisPhase = gbisParams->gbisPhase;

  //determine long and short cutoff
  float fsMax = FS_MAX; // gbisParams->fsMax;
  float a_cut = gbisParams->a_cut - fsMax;
  float a_cut_ps2 = (a_cut+fsMax)*(a_cut+fsMax);
  float r_cut = gbisParams->cutoff;
  float a_cut2 = a_cut*a_cut;
  float r_cut2 = r_cut*r_cut;
  int s = 0;// 0 is short cutoff list
  int l = 1;// 1 is long-short cutoff list
#ifdef BENCHMARK
  double t1 = 1.0*clock()/CLOCKS_PER_SEC;
  int nops = 0;
#endif
  Position ri, rj;
  Position ngri, ngrj;
  float r, dr, r2;
  float rhoi0, rhois, rhoj0, rhojs, rhois2, rhojs2, rhois2_16;
  float ami2, amj2, api2, apj2;
  int numGBISPairlists = 4;

  float max_gcut2  = r_cut;
  if (a_cut+fsMax > r_cut)
    max_gcut2 = a_cut+fsMax;
  max_gcut2 += 2.0*gbisParams->maxGroupRadius;
  //CkPrintf("max_cut = %f, %f\n",max_gcut2,gbisParams->maxGroupRadius);
  max_gcut2 *= max_gcut2;

  int maxGroupPairs = params->numAtoms[1];
  short *groupPairs = new short[maxGroupPairs];/*delete*/

  //foreach nonbonded group i
  for (int ngi = minIg; ngi < maxI; /*ngi updated at loop bottom*/ ) {
    int numGroupPairs = 0;
    ngri = params->p[0][ngi].position;
    ngri.x += offset_x;
    ngri.y += offset_y;
    ngri.z += offset_z;

    //find close j-groups; include i-group
    for (int ngj = unique*(ngi+params->p[0][ngi].nonbondedGroupSize); ngj < params->numAtoms[1]; ngj+=params->p[1][ngj].nonbondedGroupSize) {
      
      ngrj = params->p[1][ngj].position;
      dr = ngri.x - ngrj.x;
      r2 = dr*dr;
      dr = ngri.y - ngrj.y;
      r2 += dr*dr;
      dr = ngri.z - ngrj.z;
      r2 += dr*dr;
      if (r2 < max_gcut2) {
        groupPairs[numGroupPairs++] = ngj;
      }
    }

    //for all i in i-group
    int iGroupSize = params->p[0][ngi].nonbondedGroupSize;
    for (int i=ngi; i<ngi+iGroupSize; i++) {
      //CkPrintf("\tFORALL i=%05i\n",params->pExt[0][i].id);
      //extend pairlists
      int *size = new int[numGBISPairlists];
      plint **pairs = new plint*[numGBISPairlists];
      for (int k = 0; k < numGBISPairlists; k++) {
        size[k] = 0;
        pairs[k] = gbisParams->gbisStepPairlists[k]->
                             newlist(maxPairs);
      }

      //load i values
      ri = params->p[0][i].position;
      ri.x += offset_x;
      ri.y += offset_y;
      ri.z += offset_z;
      rhois = gbisParams->intRad[0][2*i+1];
      api2 = a_cut + rhois;
      api2 *= api2;
      ami2 = a_cut - rhois;
      ami2 *= ami2;
      rhois2 = rhois*rhois;
      rhois2_16 = 16.0*rhois2;

      //for all other i's in i-group
      if (unique)
      for (int j = i+1; j < ngi+iGroupSize; j++) {

#ifdef BENCHMARK
          nops ++;
#endif
          rj = params->p[1][j].position;
          dr = ri.x - rj.x;
          r2 = dr*dr;
          dr = ri.y - rj.y;
          r2 += dr*dr;
          dr = ri.z - rj.z;
          r2 += dr*dr;

          rhojs = gbisParams->intRad[1][2*j+1];
          rhojs2 = rhojs*rhojs;
          //r = sqrt(r2);
          amj2 = a_cut - rhojs;
          amj2 *= amj2;

          apj2 = a_cut + rhojs;
          apj2 *= apj2;
          //CkPrintf("IPAIR %5i %5i %5i\n",0*gbisParams->cid,params->pExt[0][i].id,params->pExt[1][j].id);

          if (  r2 < ami2 && r2 < amj2 &&
            r2 > rhois2_16 && r2 > 16.0*rhojs2 ) {
            //belongs to 22
            pairs[0][size[0]++] = j;
          } else if ( r2 < api2 && r2 < apj2 &&
            r2 > ami2 && r2 > amj2 ) {
            //belongs to 11
            pairs[1][size[1]++] = j;
          } else if ( r2 < a_cut_ps2 ) {
            //belongs to other a_cut
            pairs[2][size[2]++] = j;
          } else if ( r2 < r_cut2 ) {
            //belongs to r_cut and (r_cut > a_cut)
            pairs[3][size[3]++] = j;
          }
      }

      //foreach j-group
      for (int g = 0; g < numGroupPairs; g++) {
        int ngj = groupPairs[g];
        int jGroupSize = params->p[1][ngj].nonbondedGroupSize;
        //for each j in j-group
        int maxJ = ngj+jGroupSize;
        for (int j = ngj; j < maxJ; j++) {

          //CkPrintf("for j-atom%i\n",j);
#ifdef BENCHMARK
          nops ++;
#endif
          rj = params->p[1][j].position;
          dr = ri.x - rj.x;
          r2 = dr*dr;
          dr = ri.y - rj.y;
          r2 += dr*dr;
          dr = ri.z - rj.z;
          r2 += dr*dr;
          //CkPrintf("PAIR %5i %5i %5i\n",0*gbisParams->cid,params->pExt[0][i].id,params->pExt[1][j].id);

          rhojs = gbisParams->intRad[1][2*j+1];
          rhojs2 = rhojs*rhojs;
          //r = sqrt(r2);
          amj2 = a_cut - rhojs;
          amj2 *= amj2;

          apj2 = a_cut + rhojs;
          apj2 *= apj2;

          if (  r2 < ami2 && r2 < amj2 &&
            r2 > rhois2_16 && r2 > 16.0*rhojs2 ) {
            //belongs to 22
            //CkPrintf("Pair22 %06i %06i\n",params->pExt[0][i].id,params->pExt[1][j].id);
            pairs[0][size[0]++] = j;
          } else if ( r2 < api2 && r2 < apj2 &&
            r2 > ami2 && r2 > amj2 ) {
            //CkPrintf("Pair11 %06i %06i\n",params->pExt[0][i].id,params->pExt[1][j].id);
            //belongs to 11
            pairs[1][size[1]++] = j;
          } else if ( r2 < a_cut_ps2 ) {
            //CkPrintf("PairAA %06i %06i\n",params->pExt[0][i].id,params->pExt[1][j].id);
            //belongs to other a_cut
            pairs[2][size[2]++] = j;
          } else if ( r2 < r_cut2 ) {
            //CkPrintf("PairRR %06i %06i\n",params->pExt[0][i].id,params->pExt[1][j].id);
            //belongs to r_cut and (r_cut > a_cut)
            pairs[3][size[3]++] = j;
          }
        }//end j inner loop
      }//end j-group loop
    for (int k = 0; k < numGBISPairlists; k++) {
      gbisParams->gbisStepPairlists[k]->newsize(size[k]);
    }
    delete[] size;
    delete[] pairs;
    }//end i all atom loop

    //jump to next nbg for round-robin
    //same iteration patter followed in all three gbisPhases below
    for (int s = 0; s < strideIg; s++) {
      ngi+=params->p[0][ngi].nonbondedGroupSize;
    }
  }//end i-group loop
    delete[] groupPairs;
  for (int k = 0; k < numGBISPairlists; k++)
    gbisParams->gbisStepPairlists[k]->reset();

#ifdef BENCHMARK
  double t2 = 1.0*clock()/CLOCKS_PER_SEC;
  CkPrintf("PHASELST: %8.3f ms @ %8.3f ns/iter for %i iter\n",1000.0*(t2-t1), 1000000000.0*(t2-t1)/nops, nops);
#endif
}//end pairlist method