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,2014,2015,2016, 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 || defined LJ_EWALD)
  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, frLJ = 0;
  73             real            VLJ gmx_unused;
  74 #if defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH
  75             real            r, rsw;
  76 #endif
  77
  78 #ifdef CALC_COULOMB
  79             real qq;
  80             real fcoul;
  81 #ifdef CALC_COUL_TAB
  82             real rs, frac;
  83             int  ri;
  84             real fexcl;
  85 #endif
  86 #ifdef CALC_ENERGIES
  87             real vcoul;
  88 #endif
  89 #endif
  90             real fscal;
  91             real fx, fy, fz;
  92
  93             /* A multiply mask used to zero an interaction
  94              * when either the distance cutoff is exceeded, or
  95              * (if appropriate) the i and j indices are
  96              * unsuitable for this kind of inner loop. */
  97             real skipmask;
  98
  99 #ifdef CHECK_EXCLS
 100             /* A multiply mask used to zero an interaction
 101              * when that interaction should be excluded
 102              * (e.g. because of bonding). */
 103             int interact;
 104
 105             interact = ((l_cj[cjind].excl>>(i*UNROLLI + j)) & 1);
 106 #ifndef EXCL_FORCES
 107             skipmask = interact;
 108 #else
 109             skipmask = (cj == ci_sh && j <= i) ? 0.0 : 1.0;
 110 #endif
 111 #else
 112 #define interact 1.0
 113             skipmask = 1.0;
 114 #endif
 115
 116             // cppcheck-suppress unreadVariable
 117             VLJ = 0;
 118
 119             aj = cj*UNROLLJ + j;
 120
 121             dx  = xi[i*XI_STRIDE+XX] - x[aj*X_STRIDE+XX];
 122             dy  = xi[i*XI_STRIDE+YY] - x[aj*X_STRIDE+YY];
 123             dz  = xi[i*XI_STRIDE+ZZ] - x[aj*X_STRIDE+ZZ];
 124
 125             rsq = dx*dx + dy*dy + dz*dz;
 126
 127             /* Prepare to enforce the cut-off. */
 128             skipmask = (rsq >= rcut2) ? 0 : skipmask;
 129             /* 9 flops for r^2 + cut-off check */
 130
 131             // Ensure the distances do not fall below the limit where r^-12 overflows.
 132             // This should never happen for normal interactions.
 133             rsq = std::max(rsq, NBNXN_MIN_RSQ);
 134
 135 #ifdef COUNT_PAIRS
 136             npair++;
 137 #endif
 138
 139             rinv = gmx::invsqrt(rsq);
 140             /* 5 flops for invsqrt */
 141
 142             /* Partially enforce the cut-off (and perhaps
 143              * exclusions) to avoid possible overflow of
 144              * rinvsix when computing LJ, and/or overflowing
 145              * the Coulomb table during lookup. */
 146             rinv = rinv * skipmask;
 147
 148             rinvsq  = rinv*rinv;
 149
 150 #ifdef HALF_LJ
 151             if (i < UNROLLI/2)
 152 #endif
 153             {
 154                 c6      = nbfp[type_i_off+type[aj]*2  ];
 155                 c12     = nbfp[type_i_off+type[aj]*2+1];
 156
 157 #if defined LJ_CUT || defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH
 158                 rinvsix = interact*rinvsq*rinvsq*rinvsq;
 159                 FrLJ6   = c6*rinvsix;
 160                 FrLJ12  = c12*rinvsix*rinvsix;
 161                 frLJ    = FrLJ12 - FrLJ6;
 162                 /* 7 flops for r^-2 + LJ force */
 163 #if defined CALC_ENERGIES || defined LJ_POT_SWITCH
 164                 VLJ     = (FrLJ12 + c12*ic->repulsion_shift.cpot)/12 -
 165                     (FrLJ6 + c6*ic->dispersion_shift.cpot)/6;
 166                 /* 7 flops for LJ energy */
 167 #endif
 168 #endif
 169
 170 #if defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH
 171                 /* Force or potential switching from ic->rvdw_switch */
 172                 r       = rsq*rinv;
 173                 rsw     = r - ic->rvdw_switch;
 174                 rsw     = (rsw >= 0.0 ? rsw : 0.0);
 175 #endif
 176 #ifdef LJ_FORCE_SWITCH
 177                 frLJ   +=
 178                     -c6*(ic->dispersion_shift.c2 + ic->dispersion_shift.c3*rsw)*rsw*rsw*r
 179                     + c12*(ic->repulsion_shift.c2 + ic->repulsion_shift.c3*rsw)*rsw*rsw*r;
 180 #if defined CALC_ENERGIES
 181                 VLJ    +=
 182                     -c6*(-ic->dispersion_shift.c2/3 - ic->dispersion_shift.c3/4*rsw)*rsw*rsw*rsw
 183                     + c12*(-ic->repulsion_shift.c2/3 - ic->repulsion_shift.c3/4*rsw)*rsw*rsw*rsw;
 184 #endif
 185 #endif
 186
 187 #if defined CALC_ENERGIES || defined LJ_POT_SWITCH
 188                 /* Masking should be done after force switching,
 189                  * but before potential switching.
 190                  */
 191                 /* Need to zero the interaction if there should be exclusion. */
 192                 VLJ     = VLJ * interact;
 193 #endif
 194
 195 #ifdef LJ_POT_SWITCH
 196                 {
 197                     real sw, dsw;
 198
 199                     sw    = 1.0 + (swV3 + (swV4+ swV5*rsw)*rsw)*rsw*rsw*rsw;
 200                     dsw   = (swF2 + (swF3 + swF4*rsw)*rsw)*rsw*rsw;
 201
 202                     frLJ  = frLJ*sw - r*VLJ*dsw;
 203                     VLJ  *= sw;
 204                 }
 205 #endif
 206
 207 #ifdef LJ_EWALD
 208                 {
 209                     real            c6grid, rinvsix_nm, cr2, expmcr2, poly;
 210 #ifdef CALC_ENERGIES
 211                     real            sh_mask;
 212 #endif
 213
 214 #ifdef LJ_EWALD_COMB_GEOM
 215                     c6grid       = ljc[type[ai]*2]*ljc[type[aj]*2];
 216 #elif defined LJ_EWALD_COMB_LB
 217                     {
 218                         real sigma, sigma2, epsilon;
 219
 220                         /* These sigma and epsilon are scaled to give 6*C6 */
 221                         sigma   = ljc[type[ai]*2] + ljc[type[aj]*2];
 222                         epsilon = ljc[type[ai]*2+1]*ljc[type[aj]*2+1];
 223
 224                         sigma2  = sigma*sigma;
 225                         c6grid  = epsilon*sigma2*sigma2*sigma2;
 226                     }
 227 #else
 228 #error "No LJ Ewald combination rule defined"
 229 #endif
 230
 231 #ifdef CHECK_EXCLS
 232                     /* Recalculate rinvsix without exclusion mask */
 233                     rinvsix_nm   = rinvsq*rinvsq*rinvsq;
 234 #else
 235                     rinvsix_nm   = rinvsix;
 236 #endif
 237                     cr2          = lje_coeff2*rsq;
 238 #if GMX_DOUBLE
 239                     expmcr2      = exp(-cr2);
 240 #else
 241                     expmcr2      = expf(-cr2);
 242 #endif
 243                     poly         = 1 + cr2 + 0.5*cr2*cr2;
 244
 245                     /* Subtract the grid force from the total LJ force */
 246                     frLJ        += c6grid*(rinvsix_nm - expmcr2*(rinvsix_nm*poly + lje_coeff6_6));
 247 #ifdef CALC_ENERGIES
 248                     /* Shift should only be applied to real LJ pairs */
 249                     sh_mask      = lje_vc*interact;
 250
 251                     VLJ         += c6grid/6*(rinvsix_nm*(1 - expmcr2*poly) + sh_mask);
 252 #endif
 253                 }
 254 #endif          /* LJ_EWALD */
 255
 256 #ifdef VDW_CUTOFF_CHECK
 257                 /* Mask for VdW cut-off shorter than Coulomb cut-off */
 258                 {
 259                     real skipmask_rvdw;
 260
 261                     skipmask_rvdw = (rsq < rvdw2) ? 1.0 : 0.0;
 262                     frLJ         *= skipmask_rvdw;
 263 #ifdef CALC_ENERGIES
 264                     VLJ    *= skipmask_rvdw;
 265 #endif
 266                 }
 267 #else
 268 #if defined CALC_ENERGIES
 269                 /* Need to zero the interaction if r >= rcut */
 270                 VLJ     = VLJ * skipmask;
 271                 /* 1 more flop for LJ energy */
 272 #endif
 273 #endif          /* VDW_CUTOFF_CHECK */
 274
 275
 276 #ifdef CALC_ENERGIES
 277 #ifdef ENERGY_GROUPS
 278                 Vvdw[egp_sh_i[i]+((egp_cj>>(nbat->neg_2log*j)) & egp_mask)] += VLJ;
 279 #else
 280                 Vvdw_ci += VLJ;
 281                 /* 1 flop for LJ energy addition */
 282 #endif
 283 #endif
 284             }
 285
 286 #ifdef CALC_COULOMB
 287             /* Enforce the cut-off and perhaps exclusions. In
 288              * those cases, rinv is zero because of skipmask,
 289              * but fcoul and vcoul will later be non-zero (in
 290              * both RF and table cases) because of the
 291              * contributions that do not depend on rinv. These
 292              * contributions cannot be allowed to accumulate
 293              * to the force and potential, and the easiest way
 294              * to do this is to zero the charges in
 295              * advance. */
 296             qq = skipmask * qi[i] * q[aj];
 297
 298 #ifdef CALC_COUL_RF
 299             fcoul  = qq*(interact*rinv*rinvsq - k_rf2);
 300             /* 4 flops for RF force */
 301 #ifdef CALC_ENERGIES
 302             vcoul  = qq*(interact*rinv + k_rf*rsq - c_rf);
 303             /* 4 flops for RF energy */
 304 #endif
 305 #endif
 306
 307 #ifdef CALC_COUL_TAB
 308             rs     = rsq*rinv*ic->tabq_scale;
 309             ri     = (int)rs;
 310             frac   = rs - ri;
 311 #if !GMX_DOUBLE
 312             /* fexcl = F_i + frac * (F_(i+1)-F_i) */
 313             fexcl  = tab_coul_FDV0[ri*4] + frac*tab_coul_FDV0[ri*4+1];
 314 #else
 315             /* fexcl = (1-frac) * F_i + frac * F_(i+1) */
 316             fexcl  = (1 - frac)*tab_coul_F[ri] + frac*tab_coul_F[ri+1];
 317 #endif
 318             fcoul  = interact*rinvsq - fexcl;
 319             /* 7 flops for float 1/r-table force */
 320 #ifdef CALC_ENERGIES
 321 #if !GMX_DOUBLE
 322             vcoul  = qq*(interact*(rinv - ic->sh_ewald)
 323                          -(tab_coul_FDV0[ri*4+2]
 324                            -halfsp*frac*(tab_coul_FDV0[ri*4] + fexcl)));
 325             /* 7 flops for float 1/r-table energy (8 with excls) */
 326 #else
 327             vcoul  = qq*(interact*(rinv - ic->sh_ewald)
 328                          -(tab_coul_V[ri]
 329                            -halfsp*frac*(tab_coul_F[ri] + fexcl)));
 330 #endif
 331 #endif
 332             fcoul *= qq*rinv;
 333 #endif
 334
 335 #ifdef CALC_ENERGIES
 336 #ifdef ENERGY_GROUPS
 337             Vc[egp_sh_i[i]+((egp_cj>>(nbat->neg_2log*j)) & egp_mask)] += vcoul;
 338 #else
 339             Vc_ci += vcoul;
 340             /* 1 flop for Coulomb energy addition */
 341 #endif
 342 #endif
 343 #endif
 344
 345 #ifdef CALC_COULOMB
 346 #ifdef HALF_LJ
 347             if (i < UNROLLI/2)
 348 #endif
 349             {
 350                 fscal = frLJ*rinvsq + fcoul;
 351                 /* 2 flops for scalar LJ+Coulomb force */
 352             }
 353 #ifdef HALF_LJ
 354             else
 355             {
 356                 fscal = fcoul;
 357             }
 358 #endif
 359 #else
 360             fscal = frLJ*rinvsq;
 361 #endif
 362             fx = fscal*dx;
 363             fy = fscal*dy;
 364             fz = fscal*dz;
 365
 366             /* Increment i-atom force */
 367             fi[i*FI_STRIDE+XX] += fx;
 368             fi[i*FI_STRIDE+YY] += fy;
 369             fi[i*FI_STRIDE+ZZ] += fz;
 370             /* Decrement j-atom force */
 371             f[aj*F_STRIDE+XX]  -= fx;
 372             f[aj*F_STRIDE+YY]  -= fy;
 373             f[aj*F_STRIDE+ZZ]  -= fz;
 374             /* 9 flops for force addition */
 375         }
 376     }
 377 }
 378
 379 #undef interact
 380 #undef EXCL_FORCES