ComputeNonbondedBase2.h

Go to the documentation of this file.

EXCLUDED( FAST( foo bar ) )
EXCLUDED( MODIFIED( foo bar ) )
EXCLUDED( NORMAL( foo bar ) )
NORMAL( MODIFIED( foo bar ) )
ALCHPAIR( NOT_ALCHPAIR( foo bar ) )

ALCHPAIR(
  // get alchemical nonbonded scaling parameters (once per pairlist)
  myLambda = ALCH1(lambdaUp) ALCH2(lambdaDown) ALCH3(lambdaUp) ALCH4(lambdaDown);
  FEP(myLambda2 = ALCH1(lambda2Up) ALCH2(lambda2Down) ALCH3(lambda2Up) ALCH4(lambda2Down);)
  myElecLambda =  ALCH1(elecLambdaUp) ALCH2(elecLambdaDown) ALCH3(elecLambdaUp) ALCH4(elecLambdaDown);
  FEP(myElecLambda2 = ALCH1(elecLambda2Up) ALCH2(elecLambda2Down) ALCH3(elecLambda2Up) ALCH4(elecLambda2Down);)
  myVdwLambda =  ALCH1(vdwLambdaUp) ALCH2(vdwLambdaDown) ALCH3(vdwLambdaUp) ALCH4(vdwLambdaDown);
  FEP(myVdwLambda2 = ALCH1(vdwLambda2Up) ALCH2(vdwLambda2Down) ALCH3(vdwLambda2Up) ALCH4(vdwLambda2Down);)
  ALCH1(myRepLambda = repLambdaUp) ALCH2(myRepLambda = repLambdaDown);
  FEP(ALCH1(myRepLambda2 = repLambda2Up) ALCH2(myRepLambda2 = repLambda2Down);)
  ALCH1(myVdwShift = vdwShiftUp) ALCH2(myVdwShift = vdwShiftDown);
  FEP(ALCH1(myVdwShift2 = vdwShift2Up) ALCH2(myVdwShift2 = vdwShift2Down);)
)

#ifdef ARCH_POWERPC
     __alignx(64, table_four);
     __alignx(32, p_1);
#pragma unroll(1)
#pragma ibm independent_loop
#endif

#ifndef ARCH_POWERPC
#ifndef __INTEL_LLVM_COMPILER
#pragma ivdep
#endif
#endif


  
#if ( FULL( EXCLUDED( SHORT( 1+ ) ) ) 0 ) 
// avoid bug in Intel 15.0 compiler
#pragma novector
#else
#ifdef PRAGMA_SIMD
#ifndef TABENERGYFLAG
#ifndef GOFORCES
#pragma omp simd SHORT(FAST(reduction(+:f_i_x,f_i_y,f_i_z)) ENERGY(FAST(reduction(+:vdwEnergy) SHORT(reduction(+:electEnergy))))) \
             FULL(reduction(+:fullf_i_x,fullf_i_y,fullf_i_z) ENERGY(reduction(+:fullElectEnergy)))
#endif
#endif
#pragma loop_count avg=100
#else // PRAGMA_SIMD
#pragma loop_count avg=4
#endif // PRAGMA_SIMD
#endif
    for (k=0; k<npairi; ++k) {      
      TABENERGY(
      const int numtypes = simParams->tableNumTypes;
      const float table_spacing = simParams->tableSpacing;
      const int npertype = (int) (namdnearbyint(simParams->tableMaxDist / simParams->tableSpacing) + 1);
      )

      int table_i = (r2iilist[2*k] >> 14) + r2_delta_expc;  // table_i >= 0 
      const int j = pairlisti[k];
      //register const CompAtom *p_j = p_1 + j;
#define p_j (p_1+j)
      // const CompAtomExt *pExt_j = pExt_1 + j;
#define pExt_j_m (pExt_1+j)

      BigReal diffa = r2list[k] - r2_table[table_i];
      //const BigReal* const table_four_i = table_four + 16*table_i;
#define table_four_i (table_four + 16*table_i)

#if  ( FAST( 1 + ) TABENERGY( 1 + ) 0 ) // FAST or TABENERGY
      //const LJTable::TableEntry * lj_pars = 
      //        lj_row + 2 * p_j->vdwType MODIFIED(+ 1);
      // DJH next line we select either A,B or A14,B14 based on MODIFIED
      const int lj_index = 2 * p_j->vdwType MODIFIED(+ 1);
#define lj_pars (lj_row+lj_index)
#endif
      
      TABENERGY(
      register const int tabtype = -1 - ( lj_pars->A < 0 ? lj_pars->A : 0 );
      )
      
#if ( SHORT( FAST( 1+ ) ) 0 ) 
      //Force *f_j = f_1 + j;
#define f_j (f_1+j)
#endif
        
#if ( FULL( 1+ ) 0 )
      //Force *fullf_j = fullf_1 + j;
#define fullf_j (fullf_1+j)
#endif

      //Power PC aliasing and alignment constraints
#ifdef ARCH_POWERPC      
#if ( FULL( 1+ ) 0 )
#pragma disjoint (*table_four, *fullf_1)
#pragma disjoint (*p_1,          *fullf_1)
#pragma disjoint (*r2_table,     *fullf_1)
#pragma disjoint (*r2list,       *fullf_1)
#if ( SHORT( FAST( 1+ ) ) 0 ) 
#pragma disjoint (*f_1    ,      *fullf_1)
#pragma disjoint (*fullf_1,      *f_1)
#endif   //Short + fast
#endif   //Full

#if ( SHORT( FAST( 1+ ) ) 0 ) 
#pragma disjoint (*table_four, *f_1)
#pragma disjoint (*p_1,          *f_1)
#pragma disjoint (*r2_table,     *f_1)
#pragma disjoint (*r2list,       *f_1)
#pragma disjoint (*lj_row,      *f_1)
#endif //Short + Fast

      __alignx(64, table_four_i);
      FAST (
      __alignx(32, lj_pars);
      )
      __alignx(32, p_j);
#endif   //ARCH_POWERPC

      /*
      BigReal modf = 0.0;
      int atom2 = p_j->id;
      register char excl_flag = ( (atom2 >= excl_min && atom2 <= excl_max) ?
                                        excl_flags[atom2-excl_min] : 0 );
      if ( excl_flag ) { ++exclChecksum; }
      SELF( if ( j < j_hgroup ) { excl_flag = EXCHCK_FULL; } )
      if ( excl_flag ) {
        if ( excl_flag == EXCHCK_FULL ) {
          lj_pars = lj_null_pars;
          modf = 1.0;
        } else {
          ++lj_pars;
          modf = modf_mod;
        }
      }
      */

      BigReal kqq = kq_i * p_j->charge;
      DISP( const BigReal c6 = scaling * kd_i * pExt_j_m->dispcoef; )
      
      LES( BigReal lambda_pair = lambda_table_i[p_j->partition]; )

      register const BigReal p_ij_x = p_i_x - p_j->position.x;
      register const BigReal p_ij_y = p_i_y - p_j->position.y;
      register const BigReal p_ij_z = p_i_z - p_j->position.z;

      DISP(
        // these definitions are loop dependent
        const BigReal r2 = p_ij_x*p_ij_x + p_ij_y*p_ij_y + p_ij_z*p_ij_z;
        if (r2 >= rcut2) continue;
        const BigReal r2inv = 1/r2;
        const BigReal r6inv = r2inv * r2inv * r2inv;
        const BigReal r12inv = r6inv * r6inv;
        const BigReal aR2 = r2 * alphaLJ * alphaLJ;
        const BigReal aR4 = aR2 * aR2;
        const BigReal screen_dr = (1 - (1 + aR2 + aR4/2 + aR4*aR2/6)*exp(-aR2));
        ENERGY(
          const BigReal screen_r = (1 - (1 + aR2 + aR4/2)*exp(-aR2));
        )
      )

#if ( FAST(1+) 0 )
      const BigReal A = scaling * lj_pars->A;
      const BigReal B = scaling * lj_pars->B;
      BigReal vdw_d = A * table_four_i[0] - B * table_four_i[4];
      BigReal vdw_c = A * table_four_i[1] - B * table_four_i[5];
      BigReal vdw_b = A * table_four_i[2] - B * table_four_i[6];
      BigReal vdw_a = A * table_four_i[3] - B * table_four_i[7];

      ALCHPAIR (
        // Alchemical free energy calculation
        // Pairlist 1 and 2 are for softcore atoms, while 3 and 4 are single topology atoms. 
        // Pairlists are separated so that lambda-coupled pairs are handled
        // independently from normal nonbonded (inside ALCHPAIR macro).
        // The separation-shifted van der Waals potential and a shifted 
        // electrostatics potential for decoupling are calculated explicitly.
        // Would be faster with lookup tables but because only a small minority
        // of nonbonded pairs are lambda-coupled the impact is minimal. 
        // Explicit calculation also makes things easier to modify.
        // These are now inline functions (in ComputeNonbondedFep.C) to 
        // tidy the code
        
        const BigReal r2 = r2list[k] - r2_delta;

        // These are now inline functions (in ComputeNonbondedFep.C) to 
        // tidy the code

        FEP(
          ALCH1 (    // Don't merge/recombine the ALCH 1, 2, 3 ,4. Their functions might be modified for future algorithm changes.
          fep_vdw_forceandenergies(A, B, r2, myVdwShift, myVdwShift2,
          switchdist2, cutoff2, switchfactor, vdwForceSwitching, myVdwLambda,
          myVdwLambda2, alchWCAOn, myRepLambda, myRepLambda2, &alch_vdw_energy,
          &alch_vdw_force, &alch_vdw_energy_2);)
          ALCH2 (
          fep_vdw_forceandenergies(A, B, r2, myVdwShift, myVdwShift2,
          switchdist2, cutoff2, switchfactor, vdwForceSwitching, myVdwLambda,
          myVdwLambda2, alchWCAOn, myRepLambda, myRepLambda2, &alch_vdw_energy,
          &alch_vdw_force, &alch_vdw_energy_2);)
          ALCH3 (            // In single topology region ALCH3 & 4, all atoms are paired so softcore potential is unnecessary. 
          ENERGY(alch_vdw_energy = -myVdwLambda * (( ( diffa * vdw_d * (1/6.)+ vdw_c * (1/4.)) * diffa + vdw_b *(1/2.)) * diffa + vdw_a);)
          alch_vdw_energy_2 = -myVdwLambda2 * (( ( diffa * vdw_d * (1/6.)+ vdw_c * (1/4.)) * diffa + vdw_b *(1/2.)) * diffa + vdw_a);
          alch_vdw_force =   myVdwLambda * ((diffa * vdw_d + vdw_c) * diffa + vdw_b);)
          ALCH4 (
          ENERGY(alch_vdw_energy = -myVdwLambda * (( ( diffa * vdw_d * (1/6.)+ vdw_c * (1/4.)) * diffa + vdw_b *(1/2.)) * diffa + vdw_a);)
          alch_vdw_energy_2 = -myVdwLambda2 * (( ( diffa * vdw_d * (1/6.)+ vdw_c * (1/4.)) * diffa + vdw_b *(1/2.)) * diffa + vdw_a); 
          alch_vdw_force =   myVdwLambda * ((diffa * vdw_d + vdw_c) * diffa + vdw_b);)
          )
        TI(ti_vdw_force_energy_dUdl(A, B, r2, myVdwShift, switchdist2, cutoff2,
          switchfactor, vdwForceSwitching, myVdwLambda, alchVdwShiftCoeff,
          alchWCAOn, myRepLambda, &alch_vdw_energy, &alch_vdw_force,
          &alch_vdw_dUdl);)
      )
      
        //NOT_ALCHPAIR(
        //TABENERGY(
#if (NOT_ALCHPAIR(1+) 0)
#if (TABENERGY(1+) 0)
        if (tabtype >= 0) {
          register BigReal r1;
          r1 = sqrt(p_ij_x*p_ij_x + p_ij_y*p_ij_y + p_ij_z*p_ij_z);

          //CkPrintf("%i %i %f %f %i\n", npertype, tabtype, r1, table_spacing, (int) (namdnearbyint(r1 / table_spacing)));
          register int eneraddress;
          eneraddress = 2 * ((npertype * tabtype) + ((int) namdnearbyint(r1 / table_spacing)));
          //CkPrintf("Using distance bin %i for distance %f\n", eneraddress, r1);
          vdw_d = 0.;
          vdw_c = 0.;
          vdw_b = table_ener[eneraddress + 1] / r1;
          vdw_a = (-1/2.) * diffa * vdw_b;
          ENERGY(
            register BigReal vdw_val = table_ener[eneraddress];
            //CkPrintf("Found vdw energy of %f\n", vdw_val);
            vdwEnergy += LAM(lambda_pair *) vdw_val;
            FEP( vdwEnergy_s += d_lambda_pair * vdw_val; )
          )
        }  else {
          //)
#endif
      ENERGY(
        register BigReal vdw_val = 0;
        vdw_val -=
            ( ( diffa * vdw_d * (1/6.)+ vdw_c * (1/4.)) * diffa + vdw_b *(1/2.)) * diffa + vdw_a;
        DISP(
          vdw_val += c6 * screen_r * r6inv;
          // continuous potential at the cutoff by subtracting a constant
          BigReal potentialShift = (A*rcut12inv - B*rcut6inv);
          potentialShift += c6 * screen_rc * rcut6inv;
          vdw_val -= potentialShift;
        )

        vdwEnergy += LAM(lambda_pair *) vdw_val;

        FEP(vdwEnergy_s += vdw_val;)
      )
        //TABENERGY( } ) /* endif (tabtype >= 0) */
#if (TABENERGY (1+) 0)
        }
#endif
      //) // NOT_ALCHPAIR
#endif

      ALCHPAIR(
        ENERGY(vdwEnergy   += alch_vdw_energy;)
        FEP(vdwEnergy_s += alch_vdw_energy_2;)
        TI(ALCH1(vdwEnergy_ti_1 += alch_vdw_dUdl;) ALCH2(vdwEnergy_ti_2 += alch_vdw_dUdl;))
      ) // ALCHPAIR
      
#endif // FAST

#if ( FAST(1+) 0 )
     INT( 
      register BigReal vdw_dir;
      vdw_dir = ( diffa * vdw_d + vdw_c ) * diffa + vdw_b;
      //BigReal force_r =  LAM(lambda_pair *) vdw_dir;
      reduction[pairVDWForceIndex_X] += force_sign * vdw_dir * p_ij_x;
      reduction[pairVDWForceIndex_Y] += force_sign * vdw_dir * p_ij_y;
      reduction[pairVDWForceIndex_Z] += force_sign * vdw_dir * p_ij_z;
      )

#if ( SHORT(1+) 0 ) // Short-range electrostatics

      NORMAL(
      BigReal fast_d = kqq * table_four_i[8];
      BigReal fast_c = kqq * table_four_i[9];
      BigReal fast_b = kqq * table_four_i[10];
      BigReal fast_a = kqq * table_four_i[11];
      )
      MODIFIED(
      BigReal modfckqq = (1.0-modf_mod) * kqq;
      BigReal fast_d = modfckqq * table_four_i[8];
      BigReal fast_c = modfckqq * table_four_i[9];
      BigReal fast_b = modfckqq * table_four_i[10];
      BigReal fast_a = modfckqq * table_four_i[11];
      )

      {
      ENERGY(
      register BigReal fast_val =
        ( ( diffa * fast_d * (1/6.)+ fast_c * (1/4.)) * diffa + fast_b *(1/2.)) * diffa + fast_a;
      NOT_ALCHPAIR (
        electEnergy -=  LAM(lambda_pair *) fast_val;
        FEP(electEnergy_s -= fast_val;)
      )
      ) //ENERGY
      ALCHPAIR(
        ENERGY(electEnergy   -= myElecLambda * fast_val;)
        FEP(electEnergy_s -= myElecLambda2 * fast_val;)
        TI(
          NOENERGY(register BigReal fast_val = 
            ( ( diffa * fast_d * (1/6.)+ fast_c * (1/4.)) * diffa + fast_b *(1/2.)) * diffa + fast_a;)
          ALCH1(electEnergy_ti_1 -= fast_val;) ALCH2(electEnergy_ti_2 -= fast_val;)
        )
      )

      INT(
      register BigReal fast_dir =
      ( diffa * fast_d + fast_c ) * diffa + fast_b;
      // force_r -= -1.0 * LAM(lambda_pair *) fast_dir;
      reduction[pairElectForceIndex_X] +=  force_sign * fast_dir * p_ij_x;
      reduction[pairElectForceIndex_Y] +=  force_sign * fast_dir * p_ij_y;
      reduction[pairElectForceIndex_Z] +=  force_sign * fast_dir * p_ij_z;
      )
      }


     /*****  JE - Go  *****/
     // Now Go energy should appear in VDW place -- put vdw_b back into place
#if    ( NORMAL (1+) 0)
#if    ( GO (1+) 0)


// JLai
#ifndef CODE_REDUNDANT
#define CODE_REDUNDANT 0
#endif
#if CODE_REDUNDANT
       if (ComputeNonbondedUtil::goGroPair) {
         // Explicit goGroPair calculation; only calculates goGroPair if goGroPair is turned on
         //
         // get_gro_force has an internal checklist that sees if atom_i and atom_j are
         // in the explicit pairlist.  This is done because there is no guarantee that a 
         // processor will have atom_i and atom_j so we cannot loop over the explict atom pairs.
         // We can only loop over all pairs.
         //
         // NOTE: It does not look like fast_b is not normalized by the r vector.
         //
         // JLai
         BigReal groLJe = 0.0;
         BigReal groGausse = 0.0;
         const CompAtomExt *pExt_z = pExt_1 + j;
             BigReal groForce = mol->get_gro_force2(p_ij_x, p_ij_y, p_ij_z,pExt_i.id,pExt_z->id,&groLJe,&groGausse);
             NAMD_die("Failsafe.  This line should never be reached\n");
             fast_b += groForce;
         ENERGY(
             NOT_ALCHPAIR (
                 // JLai
                 groLJEnergy += groLJe;
                 groGaussEnergy += groGausse;
                 )
             ) //ENERGY
       }
#endif 
       BigReal goNative = 0;
       BigReal goNonnative = 0;
       BigReal goForce = 0;
       register const CompAtomExt *pExt_j = pExt_1 + j;
       if (ComputeNonbondedUtil::goMethod == 2) {
         goForce = mol->get_go_force2(p_ij_x, p_ij_y, p_ij_z, pExt_i.id, pExt_j->id,&goNative,&goNonnative);
       } else {
         //  Ported by JLai -- JE - added (
         const BigReal r2go = square(p_ij_x, p_ij_y, p_ij_z);
         const BigReal rgo = sqrt(r2go);
       
         if (ComputeNonbondedUtil::goMethod == 1) {
           goForce = mol->get_go_force(rgo, pExt_i.id, pExt_j->id, &goNative, &goNonnative);
         } else if (ComputeNonbondedUtil::goMethod == 3) {  
           goForce = mol->get_go_force_new(rgo, pExt_i.id, pExt_j->id, &goNative, &goNonnative);
         } else {
           NAMD_die("I SHOULDN'T BE HERE.  DYING MELODRAMATICALLY.\n");
         }
       }
       
       fast_b += goForce;
       {
       ENERGY(
         NOT_ALCHPAIR (
                       // JLai
                       goEnergyNative +=  goNative;
                       goEnergyNonnative += goNonnative;
         )
       ) //ENERGY                                                                                                                                                  
         INT(
           reduction[pairVDWForceIndex_X] +=  force_sign * goForce * p_ij_x;
           reduction[pairVDWForceIndex_Y] +=  force_sign * goForce * p_ij_y;
           reduction[pairVDWForceIndex_Z] +=  force_sign * goForce * p_ij_z;
         )
       }
       // End of INT 

       //DebugM(3,"rgo:" << rgo << ", pExt_i.id:" << pExt_i.id << ", pExt_j->id:" << pExt_j->id << \
         //      ", goForce:" << goForce << ", fast_b:" << fast_b << std::endl);
#endif       //     ) // End of GO macro 
       /*****  JE - End Go  *****/
       // End of port JL
#endif      //) // End of Normal MACRO

      // Combined short-range electrostatics and VdW force:
#if ( NOT_ALCHPAIR(1+) 0 )
      fast_d += vdw_d;
      fast_c += vdw_c;
      fast_b += vdw_b;
      fast_a += vdw_a;  // not used!
#endif

      register BigReal fast_dir =
                  (diffa * fast_d + fast_c) * diffa + fast_b;
      DISP(
        fast_dir += 6 * c6 * screen_dr * r6inv * r2inv;
      )
      BigReal force_r =  LAM(lambda_pair *) fast_dir;
      ALCHPAIR(
        force_r *= myElecLambda; 
        force_r += alch_vdw_force;
        // special ALCH forces already multiplied by relevant lambda
      )
          
      register BigReal tmp_x = force_r * p_ij_x;
      f_i_x += tmp_x;
      f_j->x -= tmp_x;

      register BigReal tmp_y = force_r * p_ij_y;
      f_i_y += tmp_y;
      f_j->y -= tmp_y;
      
      register BigReal tmp_z = force_r * p_ij_z;
      f_i_z += tmp_z;
      f_j->z -= tmp_z;

      PPROF(
        const BigReal p_j_z = p_j->position.z;
        int n2 = (int)floor((p_j_z-pressureProfileMin)*invThickness);
        pp_clamp(n2, pressureProfileSlabs);
        int p_j_partition = p_j->partition;

        pp_reduction(pressureProfileSlabs, n1, n2, 
                     p_i_partition, p_j_partition, pressureProfileAtomTypes,
                     tmp_x*p_ij_x, tmp_y * p_ij_y, tmp_z*p_ij_z,
                     pressureProfileReduction);
      )

#endif // SHORT
#endif // FAST

#if ( FULL (EXCLUDED( SHORT ( 1+ ) ) ) 0 ) 
      //const BigReal* const slow_i = slow_table + 4*table_i;
#define slow_i (slow_table + 4*table_i)

#ifdef ARCH_POWERPC  //Alignment and aliasing constraints
      __alignx (32, slow_table);
#if ( SHORT( FAST( 1+ ) ) 0 ) 
#pragma disjoint (*slow_table, *f_1)
#endif
#pragma disjoint (*slow_table, *fullf_1)
#endif  //ARCH_POWERPC

#endif //FULL 


#if ( FULL (MODIFIED( SHORT ( 1+ ) ) ) 0 ) 
      //const BigReal* const slow_i = slow_table + 4*table_i;
#define slow_i (slow_table + 4*table_i)

#ifdef ARCH_POWERPC //Alignment and aliasing constraints
      __alignx (32, slow_table);
#if ( SHORT( FAST( 1+ ) ) 0 ) 
#pragma disjoint (*slow_table, *f_1)
#endif
#pragma disjoint (*slow_table, *fullf_1)
#endif //ARCH_POWERPC

#endif //FULL
      
#if ( FULL( 1+ ) 0 )
      BigReal slow_d = table_four_i[8 SHORT(+ 4)];
      BigReal slow_c = table_four_i[9 SHORT(+ 4)];
      BigReal slow_b = table_four_i[10 SHORT(+ 4)];
      BigReal slow_a = table_four_i[11 SHORT(+ 4)];
      EXCLUDED(
      SHORT(
      slow_a +=    slow_i[3];
      slow_b += 2.*slow_i[2];
      slow_c += 4.*slow_i[1];
      slow_d += 6.*slow_i[0];
      )
      NOSHORT(
      slow_d -= table_four_i[12];
      slow_c -= table_four_i[13];
      slow_b -= table_four_i[14];
      slow_a -= table_four_i[15];
      )
      )
      MODIFIED(
      SHORT(
      slow_a +=    modf_mod * slow_i[3];
      slow_b += 2.*modf_mod * slow_i[2];
      slow_c += 4.*modf_mod * slow_i[1];
      slow_d += 6.*modf_mod * slow_i[0];
      )
      NOSHORT(
      slow_d -= modf_mod * table_four_i[12];
      slow_c -= modf_mod * table_four_i[13];
      slow_b -= modf_mod * table_four_i[14];
      slow_a -= modf_mod * table_four_i[15];
      )
      )
      slow_d *= kqq;
      slow_c *= kqq;
      slow_b *= kqq;
      slow_a *= kqq;

      ENERGY(
      register BigReal slow_val =
        ( ( diffa * slow_d *(1/6.)+ slow_c * (1/4.)) * diffa + slow_b *(1/2.)) * diffa + slow_a;
      
      NOT_ALCHPAIR (
        fullElectEnergy -= LAM(lambda_pair *) slow_val;
        FEP(fullElectEnergy_s -= slow_val;) 
      )
      ) // ENERGY
          
      ALCHPAIR(
        ENERGY(fullElectEnergy   -= myElecLambda * slow_val;)
        FEP(fullElectEnergy_s -= myElecLambda2 * slow_val;)
        TI(
          NOENERGY(register BigReal slow_val =
                ( ( diffa * slow_d *(1/6.)+ slow_c * (1/4.)) * diffa + slow_b *(1/2.)) * diffa + slow_a;)
          ALCH1(fullElectEnergy_ti_1 -= slow_val;) ALCH2(fullElectEnergy_ti_2 -= slow_val;)
        )
      )

      INT( {
      register BigReal slow_dir =
        ( diffa * slow_d + slow_c ) * diffa + slow_b;
      reduction[pairElectForceIndex_X] += force_sign * slow_dir * p_ij_x;
      reduction[pairElectForceIndex_Y] += force_sign * slow_dir * p_ij_y;
      reduction[pairElectForceIndex_Z] += force_sign * slow_dir * p_ij_z;
      } )


#if     (NOT_ALCHPAIR (1+) 0)
#if     (FAST(1+) 0)
#if     (NOSHORT(1+) 0)
        slow_d += vdw_d;
        slow_c += vdw_c;
        slow_b += vdw_b;
        slow_a += vdw_a; // unused!
#endif
#endif
#endif

      register BigReal slow_dir = (diffa * slow_d + slow_c) * diffa + slow_b;
      DISP(
        #if (EXCLUDED(1+) 0) || (NOSHORT(1+) 0)
        // add screened dispersion force term for either
        // excluded interaction (subtract out long-range contribution)
        // or when using merged force buffer (NOSHORT)
        slow_dir += 6 * c6 * screen_dr * r6inv * r2inv;
        #endif
        #if (EXCLUDED(1+) 0)
        ENERGY(
          // add screened dispersion energy term for
          // excluded interaction (subtract out long-range contribution)
          BigReal slow_disp_val = c6 * screen_r * r6inv;
          fullVdwEnergy += LAM(lambda_pair *) slow_disp_val;
        ) // ENERGY
        #endif
      )
      BigReal fullforce_r = slow_dir LAM(* lambda_pair);
      ALCHPAIR (
        fullforce_r *= myElecLambda;
        FAST( NOSHORT(
          fullforce_r += alch_vdw_force;
        ))
      )

      {
      register BigReal ftmp_x = fullforce_r * p_ij_x;
      fullf_i_x += ftmp_x;
      fullf_j->x -= ftmp_x;
      register BigReal ftmp_y = fullforce_r * p_ij_y;
      fullf_i_y += ftmp_y;
      fullf_j->y -= ftmp_y;
      register BigReal ftmp_z = fullforce_r * p_ij_z;
      fullf_i_z += ftmp_z;
      fullf_j->z -= ftmp_z;

      PPROF(
        const BigReal p_j_z = p_j->position.z;
        int n2 = (int)floor((p_j_z-pressureProfileMin)*invThickness);
        pp_clamp(n2, pressureProfileSlabs);
        int p_j_partition = p_j->partition;

        pp_reduction(pressureProfileSlabs, n1, n2, 
                     p_i_partition, p_j_partition, pressureProfileAtomTypes,
                     ftmp_x*p_ij_x, ftmp_y * p_ij_y, ftmp_z*p_ij_z,
                     pressureProfileReduction);

      )
      }
#endif //FULL

   } // for pairlist

#undef p_j
#undef lj_pars
#undef table_four_i
#undef slow_i
#undef f_j
#undef fullf_j