ComputeNonbondedBase.h

Go to the documentation of this file.

// Several special cases are defined:
//   NBTYPE: exclusion method (NBPAIR, NBSELF -- mutually exclusive)
//   FULLELECT full electrostatics calculation?

#ifdef ARCH_POWERPC
#include <builtins.h>
#endif

#if defined(__SSE2__) && ! defined(NAMD_DISABLE_SSE)
#include <emmintrin.h>  // We're using SSE2 intrinsics
#if defined(__INTEL_COMPILER)
#define __align(X) __declspec(align(X) )
#elif defined(__GNUC__) || defined(__PGI)
#define __align(X)  __attribute__((aligned(X) ))
#else
#define __align(X) __declspec(align(X) )
#endif
#endif

#ifdef DEFINITION // (
  #include "LJTable.h"
  #include "Molecule.h"
  #include "ComputeNonbondedUtil.h"
#endif // )

#if NAMD_ComputeNonbonded_SortAtoms != 0
  #include "PatchMap.h"
#endif

// determining class name
#undef NAME
#undef CLASS
#undef CLASSNAME
#define NAME CLASSNAME(calc)

#undef PAIR
#if NBTYPE == NBPAIR
  #define PAIR(X) X
  #define CLASS ComputeNonbondedPair
  #define CLASSNAME(X) ENERGYNAME( X ## _pair )
#else
  #define PAIR(X)
#endif

#undef SELF
#if NBTYPE == NBSELF
  #define SELF(X) X
  #define CLASS ComputeNonbondedSelf
  #define CLASSNAME(X) ENERGYNAME( X ## _self )
#else
  #define SELF(X)
#endif

#undef ENERGYNAME
#undef ENERGY
#undef NOENERGY
#ifdef CALCENERGY
  #define ENERGY(X) X
  #define NOENERGY(X)
  #define ENERGYNAME(X) SLOWONLYNAME( X ## _energy )
#else
  #define ENERGY(X)
  #define NOENERGY(X) X
  #define ENERGYNAME(X) SLOWONLYNAME( X )
#endif

#undef SLOWONLYNAME
#undef FAST
#ifdef SLOWONLY
  #define FAST(X)
  #define SLOWONLYNAME(X) MERGEELECTNAME( X ## _slow )
#else
  #define FAST(X) X
  #define SLOWONLYNAME(X) MERGEELECTNAME( X )
#endif

#undef MERGEELECTNAME
#undef SHORT
#undef NOSHORT
#ifdef MERGEELECT
  #define SHORT(X)
  #define NOSHORT(X) X
  #define MERGEELECTNAME(X) FULLELECTNAME( X ## _merge )
#else
  #define SHORT(X) X
  #define NOSHORT(X)
  #define MERGEELECTNAME(X) FULLELECTNAME( X )
#endif

#undef FULLELECTNAME
#undef FULL
#undef NOFULL
#ifdef FULLELECT
  #define FULLELECTNAME(X) TABENERGYNAME( X ## _fullelect )
  #define FULL(X) X
  #define NOFULL(X)
#else
  #define FULLELECTNAME(X) TABENERGYNAME( X )
  #define FULL(X)
  #define NOFULL(X) X
#endif

#undef TABENERGYNAME
#undef TABENERGY
#undef NOTABENERGY
#ifdef TABENERGYFLAG
  #define TABENERGYNAME(X) FEPNAME( X ## _tabener )
  #define TABENERGY(X) X
  #define NOTABENERGY(X)
#else
  #define TABENERGYNAME(X) FEPNAME( X )
  #define TABENERGY(X)
  #define NOTABENERGY(X) X
#endif

// Here are what these macros stand for:
// FEP/NOT_FEP: FEP free energy perturbation is active/inactive 
//      (does NOT use LAM)
// LES: locally-enhanced sampling is active
// LAM: scale the potential by a factor lambda (except FEP)
// INT: measure interaction energies
// PPROF: pressure profiling

#undef FEPNAME
#undef FEP
#undef LES
#undef INT
#undef PPROF
#undef LAM
#undef CUDA
#undef ALCH
#undef TI
#define FEPNAME(X) LAST( X )
#define FEP(X)
#define ALCHPAIR(X)
#define NOT_ALCHPAIR(X) X
#define LES(X)
#define INT(X)
#define PPROF(X)
#define LAM(X)
#define CUDA(X)
#define ALCH(X)
#define TI(X)
#ifdef FEPFLAG
  #undef FEPNAME
  #undef FEP
  #undef ALCH
  #define FEPNAME(X) LAST( X ## _fep )
  #define FEP(X) X
  #define ALCH(X) X
#endif
#ifdef TIFLAG
  #undef FEPNAME
  #undef TI
  #undef ALCH
  #define FEPNAME(X) LAST( X ## _ti )
  #define TI(X) X
  #define ALCH(X) X
#endif
#ifdef LESFLAG
  #undef FEPNAME
  #undef LES
  #undef LAM
  #define FEPNAME(X) LAST( X ## _les )
  #define LES(X) X
  #define LAM(X) X
#endif
#ifdef INTFLAG
  #undef FEPNAME
  #undef INT
  #define FEPNAME(X) LAST( X ## _int )
  #define INT(X) X
#endif
#ifdef PPROFFLAG
  #undef FEPNAME
  #undef INT
  #undef PPROF
  #define FEPNAME(X) LAST( X ## _pprof )
  #define INT(X) X
  #define PPROF(X) X
#endif
#ifdef NAMD_CUDA
  #undef CUDA
  #define CUDA(X) X
#endif

#define LAST(X) X

// see if things are really messed up
SELF( PAIR( foo bar ) )
LES( FEP( foo bar ) )
LES( INT( foo bar ) )
FEP( INT( foo bar ) )
LAM( INT( foo bar ) )
FEP( NOENERGY( foo bar ) )
ENERGY( NOENERGY( foo bar ) )
TABENERGY(NOTABENERGY( foo bar ) )

#if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR( + 1 ) )
  #define COMPONENT_DOTPRODUCT(A,B)  ((A##_x * B##_x) + (A##_y * B##_y) + (A##_z * B##_z))
#endif


// ************************************************************
// function header
void ComputeNonbondedUtil :: NAME
  ( nonbonded *params )

// function body
{
  // int NAME;  // to track errors via unused variable warnings

  if ( ComputeNonbondedUtil::commOnly ) return;

#ifdef NAMD_CUDA
  NOENERGY(return;)
#endif

  // speedup variables
  BigReal *reduction = params->reduction;
  SimParameters *simParams = Node::Object()->simParameters;

  PPROF(
  BigReal *pressureProfileReduction = params->pressureProfileReduction;
  const BigReal invThickness = 1.0 / pressureProfileThickness;
  )

  Pairlists &pairlists = *(params->pairlists);
#ifdef NAMD_CUDA
  int savePairlists = 0;
  int usePairlists = 0;
#else
  int savePairlists = params->savePairlists;
  int usePairlists = params->usePairlists;
#endif
  pairlists.reset();
  // PAIR(iout << "--------\n" << endi;)

  // local variables
  int exclChecksum = 0;
  FAST
  (
  ENERGY( BigReal vdwEnergy = 0; )
  SHORT
  (
  ENERGY( BigReal electEnergy = 0; )
  )

  FEP
  (
  ENERGY( BigReal vdwEnergy_s = 0; )
  SHORT
  (
  ENERGY( BigReal electEnergy_s = 0; )
  )
  )
  
  SHORT
  (
  BigReal virial_xx = 0;
  BigReal virial_xy = 0;
  BigReal virial_xz = 0;
  BigReal virial_yy = 0;
  BigReal virial_yz = 0;
  BigReal virial_zz = 0;
  )
  )
  FULL
  (
  ENERGY( BigReal fullElectEnergy = 0; )
  FEP
  (
  ENERGY( BigReal fullElectEnergy_s = 0; )
  )
  BigReal fullElectVirial_xx = 0;
  BigReal fullElectVirial_xy = 0;
  BigReal fullElectVirial_xz = 0;
  BigReal fullElectVirial_yy = 0;
  BigReal fullElectVirial_yz = 0;
  BigReal fullElectVirial_zz = 0;
  )

  // Bringing stuff into local namespace for speed.

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

  register const BigReal plcutoff2 = \
                        params->plcutoff * params->plcutoff;
  register const BigReal groupplcutoff2 = \
                        params->groupplcutoff * params->groupplcutoff;
  const BigReal dielectric_1 = ComputeNonbondedUtil:: dielectric_1;
  const LJTable* const ljTable = ComputeNonbondedUtil:: ljTable;
  LJTable::TableEntry ljNull;  ljNull.A = 0; ljNull.B = 0; TABENERGY(ljNull.tabletype=-1;)
  const LJTable::TableEntry* const lj_null_pars = &ljNull;
  const Molecule* const mol = ComputeNonbondedUtil:: mol;
  SHORT
  (
  const BigReal* const table_four = ComputeNonbondedUtil:: table_short;
  )
  FULL
  (
  SHORT
  (
  const BigReal* const slow_table = ComputeNonbondedUtil:: slow_table;
  )
  NOSHORT
  (
//#if 1 ALCH(-1)
  const BigReal* const table_four = ComputeNonbondedUtil:: table_noshort;
//#else  // have to switch this for ALCH
//  BigReal* table_four = ComputeNonbondedUtil:: table_noshort;
//#endif
  )
  )
  const BigReal scaling = ComputeNonbondedUtil:: scaling;
  const BigReal modf_mod = 1.0 - scale14;
  FAST
  (
  const BigReal switchOn2 = ComputeNonbondedUtil:: switchOn2;
  const BigReal c1 = ComputeNonbondedUtil:: c1;
  const BigReal c3 = ComputeNonbondedUtil:: c3;
  )
  const BigReal r2_delta = ComputeNonbondedUtil:: r2_delta;
  const int r2_delta_exp = ComputeNonbondedUtil:: r2_delta_exp;
  // const int r2_delta_expc = 64 * (r2_delta_exp - 127);
  const int r2_delta_expc = 64 * (r2_delta_exp - 1023);

  ALCH(
    const BigReal switchdist2 = ComputeNonbondedUtil::switchOn2;
    const BigReal cutoff2 = ComputeNonbondedUtil::cutoff2;
    const BigReal switchfactor = 1./((cutoff2 - switchdist2)*(cutoff2 - switchdist2)*(cutoff2 - switchdist2));
    const BigReal alchElecLambdaStart = ComputeNonbondedUtil::alchElecLambdaStart;
    const BigReal alchVdwLambdaEnd = ComputeNonbondedUtil::alchVdwLambdaEnd;
    const BigReal alchVdwShiftCoeff = ComputeNonbondedUtil::alchVdwShiftCoeff;

    /*lambda values 'up' are for atoms scaled up with lambda (partition 1)*/
    BigReal lambdaUp = ComputeNonbondedUtil::alchLambda;
    BigReal elecLambdaUp =  (lambdaUp <= alchElecLambdaStart)? 0. : \
              (lambdaUp - alchElecLambdaStart) / (1. - alchElecLambdaStart);
    BigReal vdwLambdaUp = 
        (lambdaUp >= alchVdwLambdaEnd)? 1. : lambdaUp / alchVdwLambdaEnd; 
    BigReal vdwShiftUp = ComputeNonbondedUtil::alchVdwShiftCoeff * (1-vdwLambdaUp);
    FEP(BigReal lambda2Up = ComputeNonbondedUtil::alchLambda2;)
    FEP(BigReal elecLambda2Up = (lambda2Up <= alchElecLambdaStart)? 0. : \
              (lambda2Up - alchElecLambdaStart) / (1. - alchElecLambdaStart);)
    FEP(BigReal vdwLambda2Up = 
        (lambda2Up >= alchVdwLambdaEnd)? 1. : lambda2Up / alchVdwLambdaEnd;) 
    FEP(BigReal vdwShift2Up = ComputeNonbondedUtil::alchVdwShiftCoeff * (1-vdwLambda2Up);)

        
    /*lambda values 'down' are for atoms scaled down with lambda (partition 2)*/
    BigReal lambdaDown = 1 - ComputeNonbondedUtil::alchLambda;
    BigReal elecLambdaDown =  (lambdaDown <= alchElecLambdaStart)? 0. : \
              (lambdaDown - alchElecLambdaStart) / (1. - alchElecLambdaStart);
    BigReal vdwLambdaDown = 
        (lambdaDown >= alchVdwLambdaEnd)? 1. : lambdaDown / alchVdwLambdaEnd; 
    BigReal vdwShiftDown = ComputeNonbondedUtil::alchVdwShiftCoeff * (1-vdwLambdaDown);
    FEP(BigReal lambda2Down = 1 - ComputeNonbondedUtil::alchLambda2;)
    FEP(BigReal elecLambda2Down = (lambda2Down <= alchElecLambdaStart)? 0. : \
        (lambda2Down - alchElecLambdaStart) / (1. - alchElecLambdaStart); )
    FEP(BigReal vdwLambda2Down = 
        (lambda2Down >= alchVdwLambdaEnd)? 1. : lambda2Down / alchVdwLambdaEnd; )
    FEP(BigReal vdwShift2Down = ComputeNonbondedUtil::alchVdwShiftCoeff * (1-vdwLambda2Down);)

  
  // Thermodynamic Integration Notes: 
  // Separation of atom pairs into different pairlist according to lambda
  // coupling, for alchemical free energy calculations. Normal pairlists
  // (pairlist[nxm]_save) are re-used for non-lambda-coupled pairs; new ones
  // (pairlist[nxm][12]_save are created for atoms switched up or down with
  // lambda respectively.
  // This makes TI coding far easier and more readable, since it's necessary 
  // to store dU/dlambda in different variables depending on whether atoms are
  // being switched up or down. Further, it allows the prior explicit coding of 
  // the separation-shifted vdW potential to stay in place without a big 
  // performance hit, since it's no longer necessary to calculate vdW potentials
  // explicity for the bulk (non-alchemical) interactions - so that part of the 
  // free energy code stays readable and easy to modify. Finally there should
  // be some performance gain over the old FEP implementation as we no
  // longer have to check the partitions of each atom pair and switch
  // parameters accordingly for every single NonbondedBase2.h loop - this is 
  // done at the pairlist level
  
  int pswitchTable[3*3];
  // determines which pairlist alchemical pairings get sent to
  // 0: non-alchemical pairs, partition 0 <-> partition 0
  // 1: atoms scaled up as lambda increases, p0<->p1
  // 2: atoms scaled down as lambda increases, p0<->p2
  // all p1<->p2 interactions to be dropped (99)
  // in general, 'up' refers to 1, 'down' refers to 2
  for (int ip=0; ip<3; ++ip) {
    for (int jp=0; jp<3; ++jp ) {
      pswitchTable[ip+3*jp] = 0;
      if ((ip==1 && jp==0) || (ip==0 && jp==1)) pswitchTable[ip+3*jp] = 1;
      if ((ip==2 && jp==0) || (ip==0 && jp==2)) pswitchTable[ip+3*jp] = 2;
      if (ip + jp == 3) pswitchTable[ip+3*jp] = 99; // no interaction

      if (! ComputeNonbondedUtil::alchDecouple) {
        // no decoupling: interactions within a partition are treated like
        // normal alchemical pairs
        if (ip == 1 && jp == 1) pswitchTable[ip+3*jp] = 1;
        if (ip == 2 && jp == 2) pswitchTable[ip+3*jp] = 2;
      }
      if (ComputeNonbondedUtil::alchDecouple) {
        // decoupling: PME calculates extra grids so that while PME 
        // interaction with the full system is switched off, a new PME grid
        // containing only alchemical atoms is switched on. Full interactions 
        // between alchemical atoms are maintained; potentials within one 
        // partition need not be scaled here.
        if (ip == 1 && jp == 1) pswitchTable[ip+3*jp] = 0;
        if (ip == 2 && jp == 2) pswitchTable[ip+3*jp] = 0;
      }
    }
  }

  // dU/dlambda variables for thermodynamic integration
  TI(
      BigReal vdwEnergy_ti_1 = 0;
      BigReal vdwEnergy_ti_2 = 0;
      SHORT(BigReal electEnergy_ti_1 = 0;
      BigReal electEnergy_ti_2 = 0;)
      FULL(BigReal fullElectEnergy_ti_1 = 0; 
      BigReal fullElectEnergy_ti_2 = 0;) 
   )
  )
        
        
  const int i_upper = params->numAtoms[0];
  register const int j_upper = params->numAtoms[1];
  register int j;
  register int i;
  const CompAtom *p_0 = params->p[0];
  const CompAtom *p_1 = params->p[1];
  const CompAtomExt *pExt_0 = params->pExt[0];
  const CompAtomExt *pExt_1 = params->pExt[1];

  char * excl_flags_buff = 0;
  const int32 * full_excl = 0;
  const int32 * mod_excl = 0;

  plint *pairlistn_save;  int npairn;
  plint *pairlistx_save;  int npairx;
  plint *pairlistm_save;  int npairm;

  ALCH(
  plint *pairlistnA1_save;  int npairnA1;
  plint *pairlistxA1_save;  int npairxA1;
  plint *pairlistmA1_save;  int npairmA1;
  plint *pairlistnA2_save;  int npairnA2;
  plint *pairlistxA2_save;  int npairxA2;
  plint *pairlistmA2_save;  int npairmA2;
  )

  NBWORKARRAYSINIT(params->workArrays);

  int arraysize = j_upper+5;

  NBWORKARRAY(plint,pairlisti,arraysize)
  NBWORKARRAY(BigReal,r2list,arraysize)

  union { double f; int32 i[2]; } byte_order_test;
  byte_order_test.f = 1.0;  // should occupy high-order bits only
  int32 *r2iilist = (int32*)r2list + ( byte_order_test.i[0] ? 0 : 1 );

  if ( ! ( savePairlists || ! usePairlists ) ) arraysize = 0;


  // DMK - Atom Sort
  // 
  // Basic Idea: Determine the line between the center of masses of
  //   the two patches.  Project and then sort the lists of atoms
  //   along this line.  Then, as the pairlists are being generated
  //   for the atoms in the first atom list, use the sorted
  //   list to only add atoms from the second list that are between
  //   +/- ~cutoff from the atoms position on the line.
  #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR( + 1 ) )

  // NOTE: For the first atom list, just the sort values themselves need to be
  // calculated (a BigReal vs. SortEntry data structure).  However, a second
  // SortEntry array is needed to perform the merge sort on the second list of
  // atoms.  Use the atomSort_[0/1]_sortValues__ arrays to sort the second
  // list of atoms and then use the left over array to make a version of the
  // list that only includes non-fixed groups (fixg).

  NBWORKARRAY(SortEntry, atomSort_0_sortValues__, arraysize);
  NBWORKARRAY(SortEntry, atomSort_1_sortValues__, arraysize);
  NBWORKARRAY(BigReal, p_0_sortValues, arraysize);

  register SortEntry* p_1_sortValues = atomSort_0_sortValues__;
  register SortEntry* p_1_sortValues_fixg = atomSort_1_sortValues__;

  int p_0_sortValues_len = 0;
  int p_1_sortValues_len = 0;
  int p_1_sortValues_fixg_len = 0;

  // Calculate the distance between to projected points on the line that
  //   represents the cutoff distance.
  BigReal atomSort_windowRadius = sqrt(groupplcutoff2);

  if (savePairlists || !usePairlists) {

    register const BigReal projLineVec_x = params->projLineVec.x;
    register const BigReal projLineVec_y = params->projLineVec.y;
    register const BigReal projLineVec_z = params->projLineVec.z;

    // Calculate the sort values for the atoms in patch 1
    {
      register int j = 0;
      register int nbgs = p_1->nonbondedGroupSize;
      register BigReal p_x = p_1->position.x;
      register BigReal p_y = p_1->position.y;
      register BigReal p_z = p_1->position.z;
      register int index = j;

      while (j < j_upper) {

        j += nbgs;

        // Set p_j_next to point to the atom for the next iteration and begin
        //   loading the 'nbgs' value for that atom.
        register const CompAtom* p_j_next = p_1 + j;
        nbgs = p_j_next->nonbondedGroupSize;

        // Calculate the distance along the projection vector
        // NOTE: If the vector from the origin to the point is 'A' and the vector
        //   between the patches is 'B' then to project 'A' onto 'B' we take the dot
        //   product of the two vectors 'A dot B' divided by the length of 'B'.
        // So... projection of A onto B = (A dot B) / length(B), however, note
        //   that length(B) == 1, therefore the sort value is simply (A dot B)
        register BigReal sortVal = COMPONENT_DOTPRODUCT(p,projLineVec);

        // Start loading the next iteration's atom's position
        p_x = p_j_next->position.x;
        p_y = p_j_next->position.y;
        p_z = p_j_next->position.z;

        // Store the caclulated sort value into the array of sort values
        register SortEntry* p_1_sortValStorePtr = p_1_sortValues + p_1_sortValues_len;
        p_1_sortValStorePtr->index = index;
        p_1_sortValStorePtr->sortValue = sortVal;
        p_1_sortValues_len++;

        // Update index for the next iteration
        index = j;       

      } // end while (j < j_upper)
    }

    // NOTE: This list and another version of it with only non-fixed
    //   atoms will be used in place of grouplist and fixglist.
    #if 0   // Selection Sort
      sortEntries_selectionSort(p_1_sortValues, p_1_sortValues_len);
    #elif 0   // Bubble Sort
      sortEntries_bubbleSort(p_1_sortValues, p_1_sortValues_len);
    #else   // Merge Sort
      #if NAMD_ComputeNonbonded_SortAtoms_LessBranches == 0
        sortEntries_mergeSort_v1(p_1_sortValues, p_1_sortValues_fixg, p_1_sortValues_len);
      #else
        sortEntries_mergeSort_v2(p_1_sortValues, p_1_sortValues_fixg, p_1_sortValues_len);
      #endif
    #endif

    // Calculate the sort values for the atoms in patch 0
    {
      register int i = 0;
      register int nbgs = p_0->nonbondedGroupSize;
      register BigReal p_x = p_0->position.x + offset_x;
      register BigReal p_y = p_0->position.y + offset_y;
      register BigReal p_z = p_0->position.z + offset_z;
      register int index = i;

      while (i < i_upper) {

        i += nbgs;

        // Set p_i_next to point to the atom for the next iteration and begin
        //   loading the 'nbgs' value for that atom.
        register const CompAtom* p_i_next = p_0 + i;
        nbgs = p_i_next->nonbondedGroupSize;

        // Calculate the distance along the projection vector
        register BigReal sortVal = COMPONENT_DOTPRODUCT(p,projLineVec);

        // Start loading the next iteration's atom's position
        p_x = p_i_next->position.x + offset_x;
        p_y = p_i_next->position.y + offset_y;
        p_z = p_i_next->position.z + offset_z;

        // Store the calculated sort value into the array of sort values
        register BigReal* p_0_sortValStorePtr = p_0_sortValues + index;
        *p_0_sortValStorePtr = sortVal;

        // Update index for the next iteration
        index = i;
      }

      p_0_sortValues_len = i_upper;
    }

  }  // end if (savePairlists || !usePairlists)

  #endif  // NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR( + 1 ) )

  // Atom Sort : The grouplist and fixglist arrays are not needed when the
  //   the atom sorting code is in use.
  #if ! (NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR( + 1 ) ) )
    NBWORKARRAY(plint,grouplist,arraysize);
    NBWORKARRAY(plint,fixglist,arraysize);
  #endif

  NBWORKARRAY(plint,goodglist,arraysize);
  NBWORKARRAY(plint,pairlistx,arraysize);
  NBWORKARRAY(plint,pairlistm,arraysize);
  NBWORKARRAY(plint,pairlist,arraysize);
  NBWORKARRAY(plint,pairlist2,arraysize);
  ALCH(
  NBWORKARRAY(plint,pairlistnAlch,arraysize);
  NBWORKARRAY(plint,pairlistnA0,arraysize);
  NBWORKARRAY(plint,pairlistnA1,arraysize);
  NBWORKARRAY(plint,pairlistxA1,arraysize);
  NBWORKARRAY(plint,pairlistmA1,arraysize);
  NBWORKARRAY(plint,pairlistnA2,arraysize);
  NBWORKARRAY(plint,pairlistxA2,arraysize);
  NBWORKARRAY(plint,pairlistmA2,arraysize);
  )

  NBWORKARRAY(short,vdwtype_array,j_upper+5);
  for (j = 0; j < j_upper; ++j)
    vdwtype_array [j] = p_1[j].vdwType;

  int fixg_upper = 0;
  int g_upper = 0;

  if ( savePairlists || ! usePairlists ) {

    #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR ( + 1 ) )

      // Create a sorted list of non-fixed groups
      register int fixg = 0;
      for (int tmpI = 0; tmpI < p_1_sortValues_len; tmpI++) {
        register SortEntry* p_1_sortEntry = p_1_sortValues + tmpI;
        register int p_1_index = p_1_sortEntry->index;
        if (!pExt_1[p_1_index].groupFixed) {
          register SortEntry* p_1_sortEntry_fixg = p_1_sortValues_fixg + p_1_sortValues_fixg_len;
          p_1_sortEntry_fixg->index = p_1_sortEntry->index;
          p_1_sortEntry_fixg->sortValue = p_1_sortEntry->sortValue;
          p_1_sortValues_fixg_len++;
        }
      }

    #else

      register int g = 0;
      for ( j = 0; j < j_upper; ++j ) {
        if ( p_1[j].nonbondedGroupSize ) {
          grouplist[g++] = j;
        }
      }
      g_upper = g;
      if ( g_upper ) grouplist[g_upper] = grouplist[g_upper-1];
      int fixg = 0;

      if ( fixedAtomsOn ) {
        for ( g = 0; g < g_upper; ++g ) {
          j = grouplist[g];
          if ( ! pExt_1[j].groupFixed ) {
            fixglist[fixg++] = j;
          }
        }
      }

      fixg_upper = fixg;
      if ( fixg_upper ) fixglist[fixg_upper] = fixglist[fixg_upper-1];

    #endif // NAMD_ComputeNonbonded_SortAtoms != 0

    *(pairlists.newlist(1)) = i_upper;
    pairlists.newsize(1);

  } else { // if ( savePairlists || ! usePairlists )

    plint *i_upper_check;
    int i_upper_check_count;
    pairlists.nextlist(&i_upper_check,&i_upper_check_count);
    if ( i_upper_check[0] != i_upper )
      NAMD_bug("pairlist i_upper mismatch!");

  } // if ( savePairlists || ! usePairlists )

  SELF(
  int j_hgroup = 0;
  int g_lower = 0;
  int fixg_lower = 0;
  )
  int pairlistindex=0;
  int pairlistoffset=0;

  

#if ( SHORT( FAST( 1+ ) ) 0 )
#if ( PAIR( 1+ ) 0 )
    Force *f_0 = params->ff[0];
    Force *f_1 = params->ff[1];
#else
#define f_1 f_0
    NBWORKARRAY(Force,f_0,i_upper)
    memset( (void*) f_0, 0, i_upper * sizeof(Force) );
#endif
#endif
#if ( FULL( 1+ ) 0 )
#if ( PAIR( 1+ ) 0 )
    Force *fullf_0 = params->fullf[0];
    Force *fullf_1 = params->fullf[1];
#else
#define fullf_1 fullf_0
    NBWORKARRAY(Force,fullf_0,i_upper);
    memset( (void*) fullf_0, 0, i_upper * sizeof(Force) );
#endif
#endif
    

  int numParts = params->numParts;
  int myPart = params->minPart;
  int groupCount = 0;

  for ( i = 0; i < (i_upper SELF(- 1)); ++i )
  {
    const CompAtom &p_i = p_0[i];
    const CompAtomExt &pExt_i = pExt_0[i];
    if ( p_i.hydrogenGroupSize ) {
      //save current group count
      int curgrpcount = groupCount;      
      //increment group count
      groupCount ++;

      if (groupCount >= numParts)
        groupCount = 0;
      
      if ( curgrpcount != myPart ) {
        i += p_i.hydrogenGroupSize - 1;
        
        //Power PC alignment constraint
#ifdef ARCH_POWERPC
        __dcbt((void *) &(p_0[i+1]));
#endif
        continue;
      }
    }

    register const BigReal p_i_x = p_i.position.x + offset_x;
    register const BigReal p_i_y = p_i.position.y + offset_y;
    register const BigReal p_i_z = p_i.position.z + offset_z;

    ALCH(const int p_i_partition = p_i.partition;)

    PPROF(
        const int p_i_partition = p_i.partition;
        int n1 = (int)floor((p_i.position.z-pressureProfileMin)*invThickness);
        pp_clamp(n1, pressureProfileSlabs);
        )

  SELF ( if ( p_i.hydrogenGroupSize ) j_hgroup = i + p_i.hydrogenGroupSize; )

  if ( savePairlists || ! usePairlists ) {

    if ( ! savePairlists ) pairlists.reset();  // limit space usage

    #ifdef MEM_OPT_VERSION
    const ExclusionCheck *exclcheck = mol->get_excl_check_for_idx(pExt_i.exclId);        
    const int excl_min = pExt_i.id + exclcheck->min;
    const int excl_max = pExt_i.id + exclcheck->max;
    #else
    const ExclusionCheck *exclcheck = mol->get_excl_check_for_atom(pExt_i.id);
    const int excl_min = exclcheck->min;
    const int excl_max = exclcheck->max;
    #endif
    const char * excl_flags_var;
    if ( exclcheck->flags ) excl_flags_var = exclcheck->flags - excl_min;
    else {  // need to build list on the fly

    //TODO: Should change later!!!!!!!!!! --Chao Mei
    //Now just for the sake of passing compilation
    #ifndef MEM_OPT_VERSION 
      if ( excl_flags_buff ) {
        int nl,l;
        nl = full_excl[0] + 1;
        for ( l=1; l<nl; ++l ) excl_flags_buff[full_excl[l]] = 0;
        nl = mod_excl[0] + 1;
        for ( l=1; l<nl; ++l ) excl_flags_buff[mod_excl[l]] = 0;
      } else {
        excl_flags_buff = new char[mol->numAtoms];
        memset( (void*) excl_flags_buff, 0, mol->numAtoms);
      }
      int nl,l;
      full_excl = mol->get_full_exclusions_for_atom(pExt_i.id);
      nl = full_excl[0] + 1;
      for ( l=1; l<nl; ++l ) excl_flags_buff[full_excl[l]] = EXCHCK_FULL;
      mod_excl = mol->get_mod_exclusions_for_atom(pExt_i.id);
      nl = mod_excl[0] + 1;
      for ( l=1; l<nl; ++l ) excl_flags_buff[mod_excl[l]] = EXCHCK_MOD;
      excl_flags_var = excl_flags_buff;
    #endif

    }
    const char * const excl_flags = excl_flags_var;

  if ( p_i.nonbondedGroupSize ) {

    pairlistindex = 0;  // initialize with 0 elements
    pairlistoffset=0;
    const int groupfixed = ( fixedAtomsOn && pExt_i.groupFixed );

    // If patch divisions are not made by hydrogen groups, then
    // hydrogenGroupSize is set to 1 for all atoms.  Thus we can
    // carry on as if we did have groups - only less efficiently.
    // An optimization in this case is to not rebuild the temporary
    // pairlist but to include every atom in it.  This should be a
    // a very minor expense.

    SELF
    (
      if ( p_i.hydrogenGroupSize ) {
        // exclude child hydrogens of i
        // j_hgroup = i + p_i.hydrogenGroupSize;  (moved above)
        while ( g_lower < g_upper &&
                grouplist[g_lower] < j_hgroup ) ++g_lower;
        while ( fixg_lower < fixg_upper &&
                fixglist[fixg_lower] < j_hgroup ) ++fixg_lower;
      }
      // add all child or sister hydrogens of i
      for ( j = i + 1; j < j_hgroup; ++j ) {
        pairlist[pairlistindex++] = j;
      }
    )

    // add remaining atoms to pairlist via hydrogen groups
    register plint *pli = pairlist + pairlistindex;

    {
      // Atom Sort : Modify the values of g and gu based on the added information
      //   of the linear projections (sort values) information.
      #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR ( + 1 ) )

        register int g = 0;
        register BigReal p_i_sortValue = p_0_sortValues[i];
        const BigReal p_i_sortValue_plus_windowRadius = p_i_sortValue + atomSort_windowRadius;
        register SortEntry* sortValues = ( groupfixed ? p_1_sortValues_fixg : p_1_sortValues );

        // Find the actual gu (upper bound in sorted list for this outer-loop atom) based on the sort values
        register int lower = 0;
        register int upper = (groupfixed ? p_1_sortValues_fixg_len : p_1_sortValues_len);
        while ((upper - lower) > 1) {
          register int j = ((lower + upper) >> 1);
          register BigReal jSortVal = sortValues[j].sortValue;
          if (jSortVal < p_i_sortValue_plus_windowRadius) {
            lower = j;
          } else {
            upper = j;
          }
        }
        const int gu = (sortValues[lower].sortValue >= p_i_sortValue_plus_windowRadius) ? lower : upper;

      #else

        register plint *gli = goodglist;
        const plint *glist = ( groupfixed ? fixglist : grouplist );
        SELF( const int gl = ( groupfixed ? fixg_lower : g_lower ); )
        const int gu = ( groupfixed ? fixg_upper : g_upper );
        register int g = PAIR(0) SELF(gl);

      #endif


      if ( g < gu ) {
        int hu = 0;
#if defined(__SSE2__) && ! defined(NAMD_DISABLE_SSE)
        if ( gu - g  >  6 ) { 

          #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR ( + 1 ) )
            register SortEntry* sortEntry0 = sortValues + g;
            register SortEntry* sortEntry1 = sortValues + g + 1;
            register int jprev0 = sortEntry0->index;
            register int jprev1 = sortEntry1->index;
          #else
            register int jprev0 = glist[g    ];
            register int jprev1 = glist[g + 1];
          #endif
          
          register int j0;
          register int j1;

          __m128d PJ_X_01 = _mm_set_pd(p_1[jprev1].position.x, p_1[jprev0].position.x);
          __m128d PJ_Y_01 = _mm_set_pd(p_1[jprev1].position.y, p_1[jprev0].position.y);
          __m128d PJ_Z_01 = _mm_set_pd(p_1[jprev1].position.z, p_1[jprev0].position.z);

          // these don't change here, so we could move them into outer scope
          const __m128d P_I_X = _mm_set1_pd(p_i_x);
          const __m128d P_I_Y = _mm_set1_pd(p_i_y);
          const __m128d P_I_Z = _mm_set1_pd(p_i_z);
 
          g += 2;
          for ( ; g < gu - 2; g +=2 ) {
            // compute 1d distance, 2-way parallel       
            j0     =  jprev0;
            j1     =  jprev1;

            __m128d T_01 = _mm_sub_pd(P_I_X, PJ_X_01);
            __m128d R2_01 = _mm_mul_pd(T_01, T_01);
            T_01 = _mm_sub_pd(P_I_Y, PJ_Y_01);
            R2_01 = _mm_add_pd(R2_01, _mm_mul_pd(T_01, T_01));
            T_01 = _mm_sub_pd(P_I_Z, PJ_Z_01);
            R2_01 = _mm_add_pd(R2_01, _mm_mul_pd(T_01, T_01));
            
            #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR ( + 1 ) )
              sortEntry0 = sortValues + g;
              sortEntry1 = sortValues + g + 1;
              jprev0 = sortEntry0->index;
              jprev1 = sortEntry1->index;
            #else
              jprev0     =  glist[g  ];
              jprev1     =  glist[g+1];
            #endif
           
            PJ_X_01 = _mm_set_pd(p_1[jprev1].position.x, p_1[jprev0].position.x);
            PJ_Y_01 = _mm_set_pd(p_1[jprev1].position.y, p_1[jprev0].position.y);
            PJ_Z_01 = _mm_set_pd(p_1[jprev1].position.z, p_1[jprev0].position.z);

            __align(16) double r2_01[2];
            _mm_store_pd(r2_01, R2_01); // 16-byte-aligned store

            // XXX these could possibly benefit from SSE-based conditionals
            bool test0 = ( r2_01[0] < groupplcutoff2 );
            bool test1 = ( r2_01[1] < groupplcutoff2 ); 
            
            //removing ifs benefits on many architectures
            //as the extra stores will only warm the cache up
            goodglist [ hu         ] = j0;
            goodglist [ hu + test0 ] = j1;
            
            hu += test0 + test1;
          }
          g-=2;
        }
#else
        if ( gu - g  >  6 ) { 

          #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR ( + 1 ) )
            register SortEntry* sortEntry0 = sortValues + g;
            register SortEntry* sortEntry1 = sortValues + g + 1;
            register int jprev0 = sortEntry0->index;
            register int jprev1 = sortEntry1->index;
          #else
            register int jprev0 = glist[g    ];
            register int jprev1 = glist[g + 1];
          #endif
          
          register  int j0; 
          register  int j1; 
          
          register  BigReal pj_x_0, pj_x_1; 
          register  BigReal pj_y_0, pj_y_1; 
          register  BigReal pj_z_0, pj_z_1; 
          register  BigReal t_0, t_1, r2_0, r2_1;
          
          pj_x_0 = p_1[jprev0].position.x;
          pj_x_1 = p_1[jprev1].position.x;  
          
          pj_y_0 = p_1[jprev0].position.y; 
          pj_y_1 = p_1[jprev1].position.y;  
          
          pj_z_0 = p_1[jprev0].position.z; 
          pj_z_1 = p_1[jprev1].position.z;
          
          g += 2;
          for ( ; g < gu - 2; g +=2 ) {
            // compute 1d distance, 2-way parallel       
            j0     =  jprev0;
            j1     =  jprev1;
            
            t_0    =  p_i_x - pj_x_0;
            t_1    =  p_i_x - pj_x_1;
            r2_0   =  t_0 * t_0;
            r2_1   =  t_1 * t_1;
            
            t_0    =  p_i_y - pj_y_0;
            t_1    =  p_i_y - pj_y_1;
            r2_0  +=  t_0 * t_0;
            r2_1  +=  t_1 * t_1;
            
            t_0    =  p_i_z - pj_z_0;
            t_1    =  p_i_z - pj_z_1;
            r2_0  +=  t_0 * t_0;
            r2_1  +=  t_1 * t_1;
            
            #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR ( + 1 ) )
              sortEntry0 = sortValues + g;
              sortEntry1 = sortValues + g + 1;
              jprev0 = sortEntry0->index;
              jprev1 = sortEntry1->index;
            #else
              jprev0     =  glist[g  ];
              jprev1     =  glist[g+1];
            #endif
            
            pj_x_0     =  p_1[jprev0].position.x;
            pj_x_1     =  p_1[jprev1].position.x;
            pj_y_0     =  p_1[jprev0].position.y; 
            pj_y_1     =  p_1[jprev1].position.y;
            pj_z_0     =  p_1[jprev0].position.z; 
            pj_z_1     =  p_1[jprev1].position.z;
            
            bool test0 = ( r2_0 < groupplcutoff2 );
            bool test1 = ( r2_1 < groupplcutoff2 ); 
            
            //removing ifs benefits on many architectures
            //as the extra stores will only warm the cache up
            goodglist [ hu         ] = j0;
            goodglist [ hu + test0 ] = j1;
            
            hu += test0 + test1;
          }
          g-=2;
        }
#endif
        
        for (; g < gu; g++) {

          #if NAMD_ComputeNonbonded_SortAtoms != 0 && ( 0 PAIR ( + 1 ) )
            register SortEntry* sortEntry = sortValues + g;
            register int j = sortEntry->index;
          #else
            int j = glist[g];
          #endif

          BigReal p_j_x = p_1[j].position.x;
          BigReal p_j_y = p_1[j].position.y;
          BigReal p_j_z = p_1[j].position.z;
          
          BigReal r2 = p_i_x - p_j_x;
          r2 *= r2;
          BigReal t2 = p_i_y - p_j_y;
          r2 += t2 * t2;
          t2 = p_i_z - p_j_z;
          r2 += t2 * t2;

          if ( r2 <= groupplcutoff2 ) 
            goodglist[hu ++] = j; 
        }

        for ( int h=0; h<hu; ++h ) {
          int j = goodglist[h];
          int nbgs = p_1[j].nonbondedGroupSize;
          pli[0] = j;   // copy over the next four in any case
          pli[1] = j+1;
          pli[2] = j+2;
          pli[3] = j+3; // assume nbgs <= 4
          pli += nbgs;
        }
      }
    }

    pairlistindex = pli - pairlist;
    // make sure padded element on pairlist points to real data
    if ( pairlistindex ) {
       pairlist[pairlistindex] = pairlist[pairlistindex-1];
    } /* PAIR( else {  // skip empty loops if no pairs were found
       int nbgs = p_i.nonbondedGroupSize;
       i += nbgs - 1;
       continue;
    } ) */
  } // if i is hydrogen group parent
  SELF
    (
    // self-comparisions require list to be incremented
    // pair-comparisions use entire list (pairlistoffset is 0)
    else pairlistoffset++;
    )

    const int atomfixed = ( fixedAtomsOn && pExt_i.atomFixed );

    register plint *pli = pairlist2;
#if 1 ALCH(-1)
    plint *pairlistn = pairlists.newlist(j_upper + 5 + 1 + 5) SELF(+ 1);
#else
    plint *pairlistn = pairlistnAlch;
#endif
    SELF( plint &pairlistn_skip = *(pairlistn-1); )
    register plint *plin = pairlistn;


    INT(
    if ( pairInteractionOn ) {
      const int ifep_type = p_i.partition;
      if (pairInteractionSelf) {
        if (ifep_type != 1) continue;
        for (int k=pairlistoffset; k<pairlistindex; k++) {
          j = pairlist[k];
          const int jfep_type = p_1[j].partition;
          // for pair-self, both atoms must be in group 1.
          if (jfep_type == 1) {
            *(pli++) = j;
          }
        }
      } else {
        if (ifep_type != 1 && ifep_type != 2) continue;
        for (int k=pairlistoffset; k<pairlistindex; k++) {
          j = pairlist[k];
          const int jfep_type = p_1[j].partition;
          // for pair, must have one from each group.
          if (ifep_type + jfep_type == 3) {
            *(pli++) = j;
          }
        }
      }
      int npair2_int = pli - pairlist2;
      pli = pairlist2;
      for (int k=0; k<npair2_int; k++) {
        j = pairlist2[k];
        BigReal p_j_x = p_1[j].position.x;
        BigReal r2 = p_i_x - p_j_x;
        r2 *= r2;
        BigReal p_j_y = p_1[j].position.y;
        BigReal t2 = p_i_y - p_j_y;
        r2 += t2 * t2;
        BigReal p_j_z = p_1[j].position.z;
        t2 = p_i_z - p_j_z;
        r2 += t2 * t2;
        if ( ( ! (atomfixed && pExt_1[j].atomFixed) ) && (r2 <= plcutoff2) ) {
          int atom2 = pExt_1[j].id;
          if ( atom2 >= excl_min && atom2 <= excl_max ) *(pli++) = j;
          else *(plin++) = j;
        }
      }
    } else
    )
    if ( atomfixed ) {
      for (int k=pairlistoffset; k<pairlistindex; k++) {
        j = pairlist[k];
        BigReal p_j_x = p_1[j].position.x;
        BigReal r2 = p_i_x - p_j_x;
        r2 *= r2;
        BigReal p_j_y = p_1[j].position.y;
        BigReal t2 = p_i_y - p_j_y;
        r2 += t2 * t2;
        BigReal p_j_z = p_1[j].position.z;
        t2 = p_i_z - p_j_z;
        r2 += t2 * t2;
        if ( (! pExt_1[j].atomFixed) && (r2 <= plcutoff2) ) {
          int atom2 = pExt_1[j].id;
          if ( atom2 >= excl_min && atom2 <= excl_max ) *(pli++) = j;
          else *(plin++) = j;
        }
      }
    } else {
      int k = pairlistoffset;
      int ku = pairlistindex;
      if ( k < ku ) {
#if defined(__SSE2__) && ! defined(NAMD_DISABLE_SSE)
        if ( ku - k  >  6 ) {      
          register  int jprev0 = pairlist [k    ];
          register  int jprev1 = pairlist [k + 1];
          
          register  int j0; 
          register  int j1; 

          __m128d PJ_X_01 = _mm_set_pd(p_1[jprev1].position.x, p_1[jprev0].position.x);
          __m128d PJ_Y_01 = _mm_set_pd(p_1[jprev1].position.y, p_1[jprev0].position.y);
          __m128d PJ_Z_01 = _mm_set_pd(p_1[jprev1].position.z, p_1[jprev0].position.z);

          // these don't change here, so we could move them into outer scope
          const __m128d P_I_X = _mm_set1_pd(p_i_x);
          const __m128d P_I_Y = _mm_set1_pd(p_i_y);
          const __m128d P_I_Z = _mm_set1_pd(p_i_z);
          
          int atom2_0 = pExt_1[jprev0].id;
          int atom2_1 = pExt_1[jprev1].id;
          
          k += 2;
          for ( ; k < ku - 2; k +=2 ) {
            // compute 1d distance, 2-way parallel       
            j0     =  jprev0;
            j1     =  jprev1;
            
            __m128d T_01 = _mm_sub_pd(P_I_X, PJ_X_01);
            __m128d R2_01 = _mm_mul_pd(T_01, T_01);
            T_01 = _mm_sub_pd(P_I_Y, PJ_Y_01);
            R2_01 = _mm_add_pd(R2_01, _mm_mul_pd(T_01, T_01));
            T_01 = _mm_sub_pd(P_I_Z, PJ_Z_01);
            R2_01 = _mm_add_pd(R2_01, _mm_mul_pd(T_01, T_01));
            
            jprev0     =  pairlist[k];
            jprev1     =  pairlist[k+1];
            
            PJ_X_01 = _mm_set_pd(p_1[jprev1].position.x, p_1[jprev0].position.x);
            PJ_Y_01 = _mm_set_pd(p_1[jprev1].position.y, p_1[jprev0].position.y);
            PJ_Z_01 = _mm_set_pd(p_1[jprev1].position.z, p_1[jprev0].position.z);

            __align(16) double r2_01[2];
            _mm_store_pd(r2_01, R2_01); // 16-byte-aligned store
            
            if (r2_01[0] <= plcutoff2) {
              if ( atom2_0 >= excl_min && atom2_0 <= excl_max ) 
                *(pli++) = j0;
              else 
                *(plin++) = j0;
            }
            atom2_0 = pExt_1[jprev0].id;
            
            if (r2_01[1] <= plcutoff2) {
              if ( atom2_1 >= excl_min && atom2_1 <= excl_max ) 
                *(pli++) = j1;
              else 
                *(plin++) = j1;
            }
            atom2_1 = pExt_1[jprev1].id;            
          }
          k-=2;
        }       
#else
        if ( ku - k  >  6 ) {      
          register  int jprev0 = pairlist [k];
          register  int jprev1 = pairlist [k + 1];
          
          register  int j0; 
          register  int j1; 
          
          register  BigReal pj_x_0, pj_x_1; 
          register  BigReal pj_y_0, pj_y_1; 
          register  BigReal pj_z_0, pj_z_1; 
          register  BigReal t_0, t_1, r2_0, r2_1;
          
          pj_x_0 = p_1[jprev0].position.x;
          pj_x_1 = p_1[jprev1].position.x;  
          
          pj_y_0 = p_1[jprev0].position.y; 
          pj_y_1 = p_1[jprev1].position.y;  
          
          pj_z_0 = p_1[jprev0].position.z; 
          pj_z_1 = p_1[jprev1].position.z;
          
          int atom2_0 = pExt_1[jprev0].id;
          int atom2_1 = pExt_1[jprev1].id;
          
          k += 2;
          for ( ; k < ku - 2; k +=2 ) {
            // compute 1d distance, 2-way parallel       
            j0     =  jprev0;
            j1     =  jprev1;
            
            t_0    =  p_i_x - pj_x_0;
            t_1    =  p_i_x - pj_x_1;
            r2_0   =  t_0 * t_0;
            r2_1   =  t_1 * t_1;
            
            t_0    =  p_i_y - pj_y_0;
            t_1    =  p_i_y - pj_y_1;
            r2_0  +=  t_0 * t_0;
            r2_1  +=  t_1 * t_1;
            
            t_0    =  p_i_z - pj_z_0;
            t_1    =  p_i_z - pj_z_1;
            r2_0  +=  t_0 * t_0;
            r2_1  +=  t_1 * t_1;
            
            jprev0     =  pairlist[k];
            jprev1     =  pairlist[k+1];
            
            pj_x_0     =  p_1[jprev0].position.x;
            pj_x_1     =  p_1[jprev1].position.x;
            pj_y_0     =  p_1[jprev0].position.y; 
            pj_y_1     =  p_1[jprev1].position.y;
            pj_z_0     =  p_1[jprev0].position.z; 
            pj_z_1     =  p_1[jprev1].position.z;
            
            if (r2_0 <= plcutoff2) {
              if ( atom2_0 >= excl_min && atom2_0 <= excl_max ) 
                *(pli++) = j0;
              else 
                *(plin++) = j0;
            }
            atom2_0 = pExt_1[jprev0].id;
            
            if (r2_1 <= plcutoff2) {
              if ( atom2_1 >= excl_min && atom2_1 <= excl_max ) 
                *(pli++) = j1;
              else 
                *(plin++) = j1;
             }
            atom2_1 = pExt_1[jprev1].id;            
          }
          k-=2;
        }       
#endif

        for (; k < ku; k++) {
          int j = pairlist[k];
          int atom2 = pExt_1[j].id;
          
          BigReal p_j_x = p_1[j].position.x;
          BigReal p_j_y = p_1[j].position.y;
          BigReal p_j_z = p_1[j].position.z;
          
          BigReal r2 = p_i_x - p_j_x;
          r2 *= r2;
          BigReal t2 = p_i_y - p_j_y;
          r2 += t2 * t2;
          t2 = p_i_z - p_j_z;
          r2 += t2 * t2;
          
          if (r2 <= plcutoff2) {
            if ( atom2 >= excl_min && atom2 <= excl_max ) 
              *(pli++) = j;
            else 
              *(plin++) = j;
          }
        }
      }
    }
    int npair2 = pli - pairlist2;
    if ( npair2 ) pairlist2[npair2] = pairlist2[npair2-1];

    plint *plix = pairlistx;
    plint *plim = pairlistm;
    plint *pln = pairlistn;
    ALCH(
    plint *plinA1 = pairlistnA1;
    plint *plixA1 = pairlistxA1;
    plint *plimA1 = pairlistmA1;
    plint *plinA2 = pairlistnA2;
    plint *plixA2 = pairlistxA2;
    plint *plimA2 = pairlistmA2;
    )
    int k=0;
#if 1 ALCH(-1)
    SELF(
    for (; pln < plin && *pln < j_hgroup; ++pln) {
      *(plix++) = *pln;  // --exclChecksum;
    }
    pairlistn_skip = pln - pairlistn;
    for (; k < npair2 && pairlist2[k] < j_hgroup; ++k) {
      *(plix++) = pairlist2[k];  // --exclChecksum;
    }
    )
#endif
ALCH(
    SELF(
    for (; pln < plin && *pln < j_hgroup; ++pln) {
      switch (pswitchTable[4*p_i_partition]) { //p_i_partition + 3*p_1[i].partition
      case 0: *(plix++) = *pln;  break;
      case 1: *(plixA1++) = *pln; break;
      case 2: *(plixA2++) = *pln; break;
      default: NAMD_die("Alchemical pairlist error"); break;
      }
    }
    for (; k < npair2 && pairlist2[k] < j_hgroup; ++k) {
      switch (pswitchTable[4*p_i_partition]) {
      case 0: *(plix++) = pairlist2[k];  break;
      case 1: *(plixA1++) = pairlist2[k]; break;
      case 2: *(plixA2++) = pairlist2[k]; break;
      default: NAMD_die("Alchemical pairlist error"); break;
      }
    }
    )
)
    for (; k < npair2; ++k ) {
      int j = pairlist2[k];
      int atom2 = pExt_1[j].id;
      int excl_flag = excl_flags[atom2];
      ALCH(int pswitch = pswitchTable[p_i_partition + 3*(p_1[j].partition)];)
      switch ( excl_flag ALCH( + 3 * pswitch)) {
      case 0:  *(plin++) = j;  break;
      case 1:  *(plix++) = j;  break;
      case 2:  *(plim++) = j;  break;
      ALCH(
      case 3:  *(plinA1++) = j; break;
      case 6:  *(plinA2++) = j; break;
      case 5:  *(plimA1++) = j; break;
      case 8:  *(plimA2++) = j; break;
      case 4:  *(plixA1++) = j; break;
      case 7:  *(plixA2++) = j; break;
      )
      }
    }
    // exclChecksum += (plix - pairlistx);
    // exclChecksum += (plim - pairlistm);


#if 1 ALCH(-1)
    npairn = plin - pln;
    pairlistn_save = pln;
    pairlistn_save[npairn] = npairn ? pairlistn_save[npairn-1] : -1;
    pairlists.newsize(plin - pairlistn SELF(+ 1) + 1);
#else
    plint *plinA0 = pairlistnA0;
    int unsortedNpairn = plin - pln;
    for ( k=0; k<unsortedNpairn; ++k ) {
      int j = pln[k];
      switch(pswitchTable[p_i_partition + 3*(p_1[j].partition)]) {
        case 0:  *(plinA0++) = j; break;
        case 1:  *(plinA1++) = j; break;
        case 2:  *(plinA2++) = j; break;
      }
    }
    
    npairn = plinA0 - pairlistnA0;
    // FB preallocation (incl extra for overhead) seems to be necessary
    pairlistn_save = pairlists.newlist(j_upper + 30);
    for ( k=0; k<npairn; ++k ) {
      pairlistn_save[k] = pairlistnA0[k];
    }
    pairlistn_save[k] = k ? pairlistn_save[k-1] : -1;
    pairlists.newsize(npairn + 1);

#endif
    
    npairx = plix - pairlistx;
    pairlistx_save = pairlists.newlist(npairx + 1);
    for ( k=0; k<npairx; ++k ) {
      pairlistx_save[k] = pairlistx[k];
    }
    pairlistx_save[k] = k ? pairlistx_save[k-1] : -1;
    pairlists.newsize(npairx + 1);

    npairm = plim - pairlistm;
    pairlistm_save = pairlists.newlist(npairm + 1);
    for ( k=0; k<npairm; ++k ) {
      pairlistm_save[k] = pairlistm[k];
    }
    pairlistm_save[k] = k ? pairlistm_save[k-1] : -1;
    pairlists.newsize(npairm + 1);

    
#if 0 ALCH(+1)
#define PAIRLISTFROMARRAY(NPAIRS,PL1,PL2,PLSAVE) \
{ \
  (NPAIRS) = (PL2) - (PL1); \
  (PLSAVE) = pairlists.newlist((NPAIRS) + 1); \
  for ( k=0; k<(NPAIRS); ++k ) { \
    (PLSAVE)[k] = (PL1)[k]; \
  } \
  (PLSAVE)[k] = k ? (PLSAVE)[k-1] : -1; \
  pairlists.newsize((NPAIRS) + 1); \
}
  
    PAIRLISTFROMARRAY(npairnA1,pairlistnA1,plinA1,pairlistnA1_save);
    PAIRLISTFROMARRAY(npairxA1,pairlistxA1,plixA1,pairlistxA1_save);
    PAIRLISTFROMARRAY(npairmA1,pairlistmA1,plimA1,pairlistmA1_save);
    PAIRLISTFROMARRAY(npairnA2,pairlistnA2,plinA2,pairlistnA2_save);
    PAIRLISTFROMARRAY(npairxA2,pairlistxA2,plixA2,pairlistxA2_save);
    PAIRLISTFROMARRAY(npairmA2,pairlistmA2,plimA2,pairlistmA2_save);
#undef PAIRLISTFROMARRAY

#endif
    
    
        // PAIR( iout << i << " " << i_upper << " save\n" << endi;)
  } else { // if ( savePairlists || ! usePairlists )
        // PAIR( iout << i << " " << i_upper << " use\n" << endi;)

    pairlists.nextlist(&pairlistn_save,&npairn);  --npairn;
    //if ( npairn > 1000 )
//      iout << i << " " << i_upper << " " << npairn << " n\n" << endi;
#if 1 ALCH(-1)
    SELF(
    int pairlistn_skip = *pairlistn_save;
    pairlistn_save += (pairlistn_skip + 1);
    npairn -= (pairlistn_skip + 1);
    )
#endif

    pairlists.nextlist(&pairlistx_save,&npairx);  --npairx;

    //if ( npairx > 1000 )
//      iout << i << " " << i_upper << " " << npairx << " x\n" << endi;
    // exclChecksum += npairx;
    pairlists.nextlist(&pairlistm_save,&npairm);  --npairm;
    ALCH(
    pairlists.nextlist(&pairlistnA1_save,&npairnA1);  --npairnA1;
    pairlists.nextlist(&pairlistxA1_save,&npairxA1);  --npairxA1;
    pairlists.nextlist(&pairlistmA1_save,&npairmA1);  --npairmA1;
    pairlists.nextlist(&pairlistnA2_save,&npairnA2);  --npairnA2;
    pairlists.nextlist(&pairlistxA2_save,&npairxA2);  --npairxA2;
    pairlists.nextlist(&pairlistmA2_save,&npairmA2);  --npairmA2;
    )
    //if ( npairm > 1000 )
//      iout << i << " " << i_upper << " " << npairm << " m\n" << endi;
    // exclChecksum += npairm;

  } // if ( savePairlists || ! usePairlists )

    LES( BigReal *lambda_table_i =
                        lambda_table + (lesFactor+1) * p_i.partition; )

    INT(
    const BigReal force_sign = (
      ( pairInteractionOn && ! pairInteractionSelf ) ?
        ( ( p_i.partition == 1 ) ? 1. : -1. ) : 0.
    );
      
    )

    const BigReal kq_i = COLOUMB * p_i.charge * scaling * dielectric_1;
    const LJTable::TableEntry * const lj_row =
                ljTable->table_row(p_i.vdwType);

    SHORT( FAST( BigReal f_i_x = 0.; ) )
    SHORT( FAST( BigReal f_i_y = 0.; ) )
    SHORT( FAST( BigReal f_i_z = 0.; ) )
    FULL( BigReal fullf_i_x = 0.; )
    FULL( BigReal fullf_i_y = 0.; )
    FULL( BigReal fullf_i_z = 0.; )

    int npairi;
    int k;

    npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
        p_i_x, p_i_y, p_i_z, p_1, pairlistn_save, npairn, pairlisti,
        r2_delta, r2list);

#define NORMAL(X) X
#define EXCLUDED(X)
#define MODIFIED(X)
#include  "ComputeNonbondedBase2.h"
#undef NORMAL
#undef EXCLUDED
#undef MODIFIED

    npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
        p_i_x, p_i_y, p_i_z, p_1, pairlistm_save, npairm, pairlisti,
        r2_delta, r2list);
    exclChecksum += npairi;

#define NORMAL(X)
#define EXCLUDED(X)
#define MODIFIED(X) X
#include  "ComputeNonbondedBase2.h"
#undef NORMAL
#undef EXCLUDED
#undef MODIFIED

#ifdef FULLELECT
    npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
        p_i_x, p_i_y, p_i_z, p_1, pairlistx_save, npairx, pairlisti,
        r2_delta, r2list);
    exclChecksum += npairi;
#ifndef NAMD_CUDA
    SELF(
    for (k=0; k<npairi && pairlisti[k] < j_hgroup; ++k) --exclChecksum;
    )
#endif

#undef FAST
#define FAST(X)
#define NORMAL(X)
#define EXCLUDED(X) X
#define MODIFIED(X)
#include  "ComputeNonbondedBase2.h"
#undef FAST
#ifdef SLOWONLY
  #define FAST(X)
#else
  #define FAST(X) X
#endif
#undef NORMAL
#undef EXCLUDED
#undef MODIFIED
#else
    exclChecksum += npairx;
#ifndef NAMD_CUDA
    SELF(
      for (k=0; k<npairx && pairlistx_save[k] < j_hgroup; ++k) --exclChecksum;
    )
#endif
#endif

#if 0 ALCH(+1)   // nonbondedbase2 for alchemical-separeted pairlists

    #undef ALCHPAIR
    #define ALCHPAIR(X) X
    #undef NOT_ALCHPAIR
    #define NOT_ALCHPAIR(X)
    BigReal myLambda; FEP(BigReal myLambda2;)
    BigReal myElecLambda;  FEP(BigReal myElecLambda2;)
    BigReal myVdwLambda; FEP(BigReal myVdwLambda2;)
    BigReal myVdwShift; FEP(BigReal myVdwShift2;)
    BigReal alch_vdw_energy; BigReal alch_vdw_force;
    FEP(BigReal alch_vdw_energy_2;) TI(BigReal alch_vdw_dUdl;)
    BigReal shiftedElec; BigReal shiftedElecForce;
    
    /********************************************************************/
    /*******NONBONDEDBASE2 FOR NORMAL INTERACTIONS SCALED BY LAMBDA******/
    /********************************************************************/
    #define NORMAL(X) X
    #define EXCLUDED(X)
    #define MODIFIED(X)

    #define ALCH1(X) X
    #define ALCH2(X)
    npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
            p_i_x, p_i_y, p_i_z, p_1, pairlistnA1_save, npairnA1, pairlisti,
            r2_delta, r2list);
    #include  "ComputeNonbondedBase2.h" // normal, direction 'up'
    #undef ALCH1
    #undef ALCH2

    #define ALCH1(X)
    #define ALCH2(X) X
    npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
            p_i_x, p_i_y, p_i_z, p_1, pairlistnA2_save, npairnA2, pairlisti,
            r2_delta, r2list);
    #include  "ComputeNonbondedBase2.h" // normal, direction 'down'
    #undef ALCH1
    #undef ALCH2

    #undef NORMAL
    #undef EXCLUDED
    #undef MODIFIED
    
    /********************************************************************/
    /******NONBONDEDBASE2 FOR MODIFIED INTERACTIONS SCALED BY LAMBDA*****/
    /********************************************************************/
    #define NORMAL(X)
    #define EXCLUDED(X)
    #define MODIFIED(X) X

    #define ALCH1(X) X
    #define ALCH2(X)
    npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
            p_i_x, p_i_y, p_i_z, p_1, pairlistmA1_save, npairmA1, pairlisti,
            r2_delta, r2list);
        exclChecksum += npairi;
    #include  "ComputeNonbondedBase2.h" // modified, direction 'up'
    #undef ALCH1
    #undef ALCH2

    #define ALCH1(X)
    #define ALCH2(X) X
    npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
            p_i_x, p_i_y, p_i_z, p_1, pairlistmA2_save, npairmA2, pairlisti,
            r2_delta, r2list);
        exclChecksum += npairi;
    #include  "ComputeNonbondedBase2.h" // modified, direction 'down'
    #undef ALCH1
    #undef ALCH2

    #undef NORMAL
    #undef EXCLUDED
    #undef MODIFIED
    
    /********************************************************************/
    /******NONBONDEDBASE2 FOR EXCLUDED INTERACTIONS SCALED BY LAMBDA*****/
    /********************************************************************/
    #ifdef FULLELECT
    #undef FAST
    #define FAST(X)
    #define NORMAL(X)
    #define EXCLUDED(X) X
    #define MODIFIED(X)

    #define ALCH1(X) X
    #define ALCH2(X)
    npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
            p_i_x, p_i_y, p_i_z, p_1, pairlistxA1_save, npairxA1, pairlisti,
            r2_delta, r2list);
        exclChecksum += npairi;
        SELF(
        for (k=0; k<npairi && pairlisti[k] < j_hgroup; ++k) --exclChecksum;
        )
    #include  "ComputeNonbondedBase2.h"  //excluded, direction 'up'
    #undef ALCH1
    #undef ALCH2

    #define ALCH1(X)
    #define ALCH2(X) X
    npairi = pairlist_from_pairlist(ComputeNonbondedUtil::cutoff2,
            p_i_x, p_i_y, p_i_z, p_1, pairlistxA2_save, npairxA2, pairlisti,
            r2_delta, r2list);
        exclChecksum += npairi;
        SELF(
        for (k=0; k<npairi && pairlisti[k] < j_hgroup; ++k) --exclChecksum;
        )
    #include  "ComputeNonbondedBase2.h"  //excluded, direction 'down'
    #undef ALCH1
    #undef ALCH2

    #undef FAST
    #ifdef SLOWONLY
      #define FAST(X)
    #else
      #define FAST(X) X
    #endif
    #undef NORMAL
    #undef EXCLUDED
    #undef MODIFIED
    #else
        exclChecksum += npairxA1 + npairxA2;
        SELF(
          for (k=0; k<npairxA1 && pairlistxA1_save[k] < j_hgroup; ++k) --exclChecksum;
          for (k=0; k<npairxA2 && pairlistxA2_save[k] < j_hgroup; ++k) --exclChecksum;
        )
    #endif

    #undef ALCHPAIR
    #define ALCHPAIR(X)
    #undef NOT_ALCHPAIR
    #define NOT_ALCHPAIR(X) X

#endif // end nonbondedbase2.h includes for alchemical pairlists

#ifdef NETWORK_PROGRESS
    CkNetworkProgress();
#endif

#ifdef ARCH_POWERPC
    //data cache block touch the position structure
    __dcbt ((void *) &(p_0[i+1]));
    __prefetch_by_load ((void *)&(groupCount));
#endif

#ifndef NAMD_CUDA
    SHORT( FAST( f_0[i].x += f_i_x; ) )
    SHORT( FAST( f_0[i].y += f_i_y; ) )
    SHORT( FAST( f_0[i].z += f_i_z; ) )
    FULL( fullf_0[i].x += fullf_i_x; )
    FULL( fullf_0[i].y += fullf_i_y; )
    FULL( fullf_0[i].z += fullf_i_z; )
#endif

        // PAIR( iout << i << " " << i_upper << " end\n" << endi;)
  } // for i

  // PAIR(iout << "++++++++\n" << endi;)

#ifdef f_1
#undef f_1
#endif
#if ( SHORT( FAST( 1+ ) ) 0 )
#if ( SELF( 1+ ) 0 )
  {
    Force *patch_f_0 = params->ff[0];

#ifndef NAMD_CUDA
#ifndef ARCH_POWERPC 
#pragma ivdep
#endif
    for ( int i = 0; i < i_upper; ++i ) {
      patch_f_0[i].x += f_0[i].x;
      patch_f_0[i].y += f_0[i].y;
      patch_f_0[i].z += f_0[i].z;
      virial_xx += f_0[i].x * p_0[i].position.x;
      virial_xy += f_0[i].x * p_0[i].position.y;
      virial_xz += f_0[i].x * p_0[i].position.z;
      virial_yy += f_0[i].y * p_0[i].position.y;
      virial_yz += f_0[i].y * p_0[i].position.z;
      virial_zz += f_0[i].z * p_0[i].position.z;
    }
#endif
  }
#endif
#endif
#ifdef fullf_1
#undef fullf_1
#endif
#if ( FULL( 1+ ) 0 )
#if ( SELF( 1+ ) 0 )
  {
    Force *patch_fullf_0 = params->fullf[0];
#ifndef NAMD_CUDA
#ifndef ARCH_POWERPC 
#pragma ivdep
#endif
    for ( int i = 0; i < i_upper; ++i ) {
      patch_fullf_0[i].x += fullf_0[i].x;
      patch_fullf_0[i].y += fullf_0[i].y;
      patch_fullf_0[i].z += fullf_0[i].z;
      fullElectVirial_xx += fullf_0[i].x * p_0[i].position.x;
      fullElectVirial_xy += fullf_0[i].x * p_0[i].position.y;
      fullElectVirial_xz += fullf_0[i].x * p_0[i].position.z;
      fullElectVirial_yy += fullf_0[i].y * p_0[i].position.y;
      fullElectVirial_yz += fullf_0[i].y * p_0[i].position.z;
      fullElectVirial_zz += fullf_0[i].z * p_0[i].position.z;
    }
#endif
  }
#endif
#endif

#ifndef NAMD_CUDA
  reduction[exclChecksumIndex] += exclChecksum;
#endif
  FAST
  (
  ENERGY( reduction[vdwEnergyIndex] += vdwEnergy; )
  SHORT
  (
  ENERGY( reduction[electEnergyIndex] += electEnergy; )
  )
  ALCH
  (
    TI(reduction[vdwEnergyIndex_ti_1] += vdwEnergy_ti_1;) 
    TI(reduction[vdwEnergyIndex_ti_2] += vdwEnergy_ti_2;) 
    FEP( reduction[vdwEnergyIndex_s] += vdwEnergy_s; )
  SHORT
  (
    FEP( reduction[electEnergyIndex_s] += electEnergy_s; )
    TI(reduction[electEnergyIndex_ti_1] += electEnergy_ti_1;) 
    TI(reduction[electEnergyIndex_ti_2] += electEnergy_ti_2;) 
  )
  )
  SHORT
  (
  reduction[virialIndex_XX] += virial_xx;
  reduction[virialIndex_XY] += virial_xy;
  reduction[virialIndex_XZ] += virial_xz;
  reduction[virialIndex_YX] += virial_xy;
  reduction[virialIndex_YY] += virial_yy;
  reduction[virialIndex_YZ] += virial_yz;
  reduction[virialIndex_ZX] += virial_xz;
  reduction[virialIndex_ZY] += virial_yz;
  reduction[virialIndex_ZZ] += virial_zz;
  )
  )
  FULL
  (
  ENERGY( reduction[fullElectEnergyIndex] += fullElectEnergy; )
  ALCH
  (
    FEP( reduction[fullElectEnergyIndex_s] += fullElectEnergy_s; )
    TI(reduction[fullElectEnergyIndex_ti_1] += fullElectEnergy_ti_1;) 
    TI(reduction[fullElectEnergyIndex_ti_2] += fullElectEnergy_ti_2;) 
  )
  reduction[fullElectVirialIndex_XX] += fullElectVirial_xx;
  reduction[fullElectVirialIndex_XY] += fullElectVirial_xy;
  reduction[fullElectVirialIndex_XZ] += fullElectVirial_xz;
  reduction[fullElectVirialIndex_YX] += fullElectVirial_xy;
  reduction[fullElectVirialIndex_YY] += fullElectVirial_yy;
  reduction[fullElectVirialIndex_YZ] += fullElectVirial_yz;
  reduction[fullElectVirialIndex_ZX] += fullElectVirial_xz;
  reduction[fullElectVirialIndex_ZY] += fullElectVirial_yz;
  reduction[fullElectVirialIndex_ZZ] += fullElectVirial_zz;
  )

  delete [] excl_flags_buff;

}

---