src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref_inner.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 /* When calculating RF or Ewald interactions we calculate the electrostatic
  37  * forces and energies on excluded atom pairs here in the non-bonded loops.
  38  */
  39 #if defined CHECK_EXCLS && defined CALC_COULOMB
  40 #define EXCL_FORCES
  41 #endif
  42
  43 {
  44     int cj;
  45 #ifdef ENERGY_GROUPS
  46     int egp_cj;
  47 #endif
  48     int i;
  49
  50     cj = l_cj[cjind].cj;
  51
  52 #ifdef ENERGY_GROUPS
  53     egp_cj = nbat->energrp[cj];
  54 #endif
  55     for (i = 0; i < UNROLLI; i++)
  56     {
  57         int ai;
  58         int type_i_off;
  59         int j;
  60
  61         ai = ci*UNROLLI + i;
  62
  63         type_i_off = type[ai]*ntype2;
  64
  65         for (j = 0; j < UNROLLJ; j++)
  66         {
  67             int  aj;
  68             real dx, dy, dz;
  69             real rsq, rinv;
  70             real rinvsq, rinvsix;
  71             real c6, c12;
  72             real FrLJ6 = 0, FrLJ12 = 0, VLJ = 0;
  73 #ifdef CALC_COULOMB
  74             real qq;
  75             real fcoul;
  76 #ifdef CALC_COUL_TAB
  77             real rs, frac;
  78             int  ri;
  79             real fexcl;
  80 #endif
  81 #ifdef CALC_ENERGIES
  82             real vcoul;
  83 #endif
  84 #endif
  85             real fscal;
  86             real fx, fy, fz;
  87
  88             /* A multiply mask used to zero an interaction
  89              * when either the distance cutoff is exceeded, or
  90              * (if appropriate) the i and j indices are
  91              * unsuitable for this kind of inner loop. */
  92             real skipmask;
  93 #ifdef VDW_CUTOFF_CHECK
  94             real skipmask_rvdw;
  95 #endif
  96 #ifdef CHECK_EXCLS
  97             /* A multiply mask used to zero an interaction
  98              * when that interaction should be excluded
  99              * (e.g. because of bonding). */
 100             int interact;
 101
 102             interact = ((l_cj[cjind].excl>>(i*UNROLLI + j)) & 1);
 103 #ifndef EXCL_FORCES
 104             skipmask = interact;
 105 #else
 106             skipmask = !(cj == ci_sh && j <= i);
 107 #endif
 108 #else
 109 #define interact 1.0
 110             skipmask = 1.0;
 111 #endif
 112
 113             aj = cj*UNROLLJ + j;
 114
 115             dx  = xi[i*XI_STRIDE+XX] - x[aj*X_STRIDE+XX];
 116             dy  = xi[i*XI_STRIDE+YY] - x[aj*X_STRIDE+YY];
 117             dz  = xi[i*XI_STRIDE+ZZ] - x[aj*X_STRIDE+ZZ];
 118
 119             rsq = dx*dx + dy*dy + dz*dz;
 120
 121             /* Prepare to enforce the cut-off. */
 122             skipmask = (rsq >= rcut2) ? 0 : skipmask;
 123             /* 9 flops for r^2 + cut-off check */
 124
 125 #ifdef CHECK_EXCLS
 126             /* Excluded atoms are allowed to be on top of each other.
 127              * To avoid overflow of rinv, rinvsq and rinvsix
 128              * we add a small number to rsq for excluded pairs only.
 129              */
 130             rsq += (1 - interact)*NBNXN_AVOID_SING_R2_INC;
 131 #endif
 132
 133 #ifdef COUNT_PAIRS
 134             npair++;
 135 #endif
 136
 137             rinv = gmx_invsqrt(rsq);
 138             /* 5 flops for invsqrt */
 139
 140             /* Partially enforce the cut-off (and perhaps
 141              * exclusions) to avoid possible overflow of
 142              * rinvsix when computing LJ, and/or overflowing
 143              * the Coulomb table during lookup. */
 144             rinv = rinv * skipmask;
 145
 146             rinvsq  = rinv*rinv;
 147
 148 #ifdef HALF_LJ
 149             if (i < UNROLLI/2)
 150 #endif
 151             {
 152                 rinvsix = interact*rinvsq*rinvsq*rinvsq;
 153
 154 #ifdef VDW_CUTOFF_CHECK
 155                 skipmask_rvdw = (rsq < rvdw2);
 156                 rinvsix      *= skipmask_rvdw;
 157 #endif
 158
 159                 c6      = nbfp[type_i_off+type[aj]*2  ];
 160                 c12     = nbfp[type_i_off+type[aj]*2+1];
 161                 FrLJ6   = c6*rinvsix;
 162                 FrLJ12  = c12*rinvsix*rinvsix;
 163                 /* 6 flops for r^-2 + LJ force */
 164 #ifdef CALC_ENERGIES
 165                 VLJ     = (FrLJ12 - c12*sh_invrc6*sh_invrc6)/12 -
 166                     (FrLJ6 - c6*sh_invrc6)/6;
 167                 /* Need to zero the interaction if r >= rcut
 168                  * or there should be exclusion. */
 169                 VLJ     = VLJ * skipmask * interact;
 170                 /* 9 flops for LJ energy */
 171 #ifdef VDW_CUTOFF_CHECK
 172                 VLJ    *= skipmask_rvdw;
 173 #endif
 174 #ifdef ENERGY_GROUPS
 175                 Vvdw[egp_sh_i[i]+((egp_cj>>(nbat->neg_2log*j)) & egp_mask)] += VLJ;
 176 #else
 177                 Vvdw_ci += VLJ;
 178                 /* 1 flop for LJ energy addition */
 179 #endif
 180 #endif
 181             }
 182
 183 #ifdef CALC_COULOMB
 184             /* Enforce the cut-off and perhaps exclusions. In
 185              * those cases, rinv is zero because of skipmask,
 186              * but fcoul and vcoul will later be non-zero (in
 187              * both RF and table cases) because of the
 188              * contributions that do not depend on rinv. These
 189              * contributions cannot be allowed to accumulate
 190              * to the force and potential, and the easiest way
 191              * to do this is to zero the charges in
 192              * advance. */
 193             qq = skipmask * qi[i] * q[aj];
 194
 195 #ifdef CALC_COUL_RF
 196             fcoul  = qq*(interact*rinv*rinvsq - k_rf2);
 197             /* 4 flops for RF force */
 198 #ifdef CALC_ENERGIES
 199             vcoul  = qq*(interact*rinv + k_rf*rsq - c_rf);
 200             /* 4 flops for RF energy */
 201 #endif
 202 #endif
 203
 204 #ifdef CALC_COUL_TAB
 205             rs     = rsq*rinv*ic->tabq_scale;
 206             ri     = (int)rs;
 207             frac   = rs - ri;
 208 #ifndef GMX_DOUBLE
 209             /* fexcl = F_i + frac * (F_(i+1)-F_i) */
 210             fexcl  = tab_coul_FDV0[ri*4] + frac*tab_coul_FDV0[ri*4+1];
 211 #else
 212             /* fexcl = (1-frac) * F_i + frac * F_(i+1) */
 213             fexcl  = (1 - frac)*tab_coul_F[ri] + frac*tab_coul_F[ri+1];
 214 #endif
 215             fcoul  = interact*rinvsq - fexcl;
 216             /* 7 flops for float 1/r-table force */
 217 #ifdef CALC_ENERGIES
 218 #ifndef GMX_DOUBLE
 219             vcoul  = qq*(interact*(rinv - ic->sh_ewald)
 220                          -(tab_coul_FDV0[ri*4+2]
 221                            -halfsp*frac*(tab_coul_FDV0[ri*4] + fexcl)));
 222             /* 7 flops for float 1/r-table energy (8 with excls) */
 223 #else
 224             vcoul  = qq*(interact*(rinv - ic->sh_ewald)
 225                          -(tab_coul_V[ri]
 226                            -halfsp*frac*(tab_coul_F[ri] + fexcl)));
 227 #endif
 228 #endif
 229             fcoul *= qq*rinv;
 230 #endif
 231
 232 #ifdef CALC_ENERGIES
 233 #ifdef ENERGY_GROUPS
 234             Vc[egp_sh_i[i]+((egp_cj>>(nbat->neg_2log*j)) & egp_mask)] += vcoul;
 235 #else
 236             Vc_ci += vcoul;
 237             /* 1 flop for Coulomb energy addition */
 238 #endif
 239 #endif
 240 #endif
 241
 242 #ifdef CALC_COULOMB
 243 #ifdef HALF_LJ
 244             if (i < UNROLLI/2)
 245 #endif
 246             {
 247                 fscal = (FrLJ12 - FrLJ6)*rinvsq + fcoul;
 248                 /* 3 flops for scalar LJ+Coulomb force */
 249             }
 250 #ifdef HALF_LJ
 251             else
 252             {
 253                 fscal = fcoul;
 254             }
 255 #endif
 256 #else
 257             fscal = (FrLJ12 - FrLJ6)*rinvsq;
 258 #endif
 259             fx = fscal*dx;
 260             fy = fscal*dy;
 261             fz = fscal*dz;
 262
 263             /* Increment i-atom force */
 264             fi[i*FI_STRIDE+XX] += fx;
 265             fi[i*FI_STRIDE+YY] += fy;
 266             fi[i*FI_STRIDE+ZZ] += fz;
 267             /* Decrement j-atom force */
 268             f[aj*F_STRIDE+XX]  -= fx;
 269             f[aj*F_STRIDE+YY]  -= fy;
 270             f[aj*F_STRIDE+ZZ]  -= fz;
 271             /* 9 flops for force addition */
 272         }
 273     }
 274 }
 275
 276 #undef interact
 277 #undef EXCL_FORCES