src/gromacs/nbnxm/kernels_reference/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.
   5  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36
  37 /* When calculating RF or Ewald interactions we calculate the electrostatic
  38  * forces and energies on excluded atom pairs here in the non-bonded loops.
  39  */
  40 #if defined CHECK_EXCLS && (defined CALC_COULOMB || defined LJ_EWALD)
  41 #    define EXCL_FORCES
  42 #endif
  43
  44 {
  45     int cj;
  46 #ifdef ENERGY_GROUPS
  47     int egp_cj;
  48 #endif
  49     int i;
  50
  51     cj = l_cj[cjind].cj;
  52
  53 #ifdef ENERGY_GROUPS
  54     egp_cj = nbatParams.energrp[cj];
  55 #endif
  56     for (i = 0; i < UNROLLI; i++)
  57     {
  58         int ai;
  59         int type_i_off;
  60         int j;
  61
  62         ai = ci * UNROLLI + i;
  63
  64         type_i_off = type[ai] * ntype2;
  65
  66         for (j = 0; j < UNROLLJ; j++)
  67         {
  68             int  aj;
  69             real dx, dy, dz;
  70             real rsq, rinv;
  71             real rinvsq, rinvsix;
  72             real c6, c12;
  73             real FrLJ6 = 0, FrLJ12 = 0, frLJ = 0;
  74             real VLJ gmx_unused;
  75 #if defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH
  76             real r, rsw;
  77 #endif
  78
  79 #ifdef CALC_COULOMB
  80             real qq;
  81             real fcoul;
  82 #    ifdef CALC_COUL_TAB
  83             real rs, frac;
  84             int  ri;
  85             real fexcl;
  86 #    endif
  87 #    ifdef CALC_ENERGIES
  88             real vcoul;
  89 #    endif
  90 #endif
  91             real fscal;
  92             real fx, fy, fz;
  93
  94             /* A multiply mask used to zero an interaction
  95              * when either the distance cutoff is exceeded, or
  96              * (if appropriate) the i and j indices are
  97              * unsuitable for this kind of inner loop. */
  98             real skipmask;
  99
 100 #ifdef CHECK_EXCLS
 101             /* A multiply mask used to zero an interaction
 102              * when that interaction should be excluded
 103              * (e.g. because of bonding). */
 104             real interact = static_cast<real>((l_cj[cjind].excl >> (i * UNROLLI + j)) & 1);
 105 #    ifndef EXCL_FORCES
 106             skipmask = interact;
 107 #    else
 108             skipmask = (cj == ci_sh && j <= i) ? 0.0 : 1.0;
 109 #    endif
 110 #else
 111             constexpr real interact = 1.0;
 112             skipmask                = interact;
 113 #endif
 114
 115             VLJ = 0;
 116
 117             aj = cj * UNROLLJ + j;
 118
 119             dx = xi[i * XI_STRIDE + XX] - x[aj * X_STRIDE + XX];
 120             dy = xi[i * XI_STRIDE + YY] - x[aj * X_STRIDE + YY];
 121             dz = xi[i * XI_STRIDE + ZZ] - x[aj * X_STRIDE + ZZ];
 122
 123             rsq = dx * dx + dy * dy + dz * dz;
 124
 125             /* Prepare to enforce the cut-off. */
 126             skipmask = (rsq >= rcut2) ? 0 : skipmask;
 127             /* 9 flops for r^2 + cut-off check */
 128
 129             // Ensure the distances do not fall below the limit where r^-12 overflows.
 130             // This should never happen for normal interactions.
 131             rsq = std::max(rsq, c_nbnxnMinDistanceSquared);
 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                 c6  = nbfp[type_i_off + type[aj] * 2];
 153                 c12 = nbfp[type_i_off + type[aj] * 2 + 1];
 154
 155 #if defined LJ_CUT || defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH
 156                 rinvsix = interact * rinvsq * rinvsq * rinvsq;
 157                 FrLJ6   = c6 * rinvsix;
 158                 FrLJ12  = c12 * rinvsix * rinvsix;
 159                 frLJ    = FrLJ12 - FrLJ6;
 160                 /* 7 flops for r^-2 + LJ force */
 161 #    if defined CALC_ENERGIES || defined LJ_POT_SWITCH
 162                 VLJ = (FrLJ12 + c12 * ic->repulsion_shift.cpot) / 12
 163                       - (FrLJ6 + c6 * ic->dispersion_shift.cpot) / 6;
 164                 /* 7 flops for LJ energy */
 165 #    endif
 166 #endif
 167
 168 #if defined LJ_FORCE_SWITCH || defined LJ_POT_SWITCH
 169                 /* Force or potential switching from ic->rvdw_switch */
 170                 r   = rsq * rinv;
 171                 rsw = r - ic->rvdw_switch;
 172                 rsw = (rsw >= 0.0 ? rsw : 0.0);
 173 #endif
 174 #ifdef LJ_FORCE_SWITCH
 175                 frLJ += -c6 * (ic->dispersion_shift.c2 + ic->dispersion_shift.c3 * rsw) * rsw * rsw * r
 176                         + c12 * (ic->repulsion_shift.c2 + ic->repulsion_shift.c3 * rsw) * rsw * rsw * r;
 177 #    if defined CALC_ENERGIES
 178                 VLJ += -c6 * (-ic->dispersion_shift.c2 / 3 - ic->dispersion_shift.c3 / 4 * rsw)
 179                                * rsw * rsw * rsw
 180                        + c12 * (-ic->repulsion_shift.c2 / 3 - ic->repulsion_shift.c3 / 4 * rsw)
 181                                  * rsw * rsw * rsw;
 182 #    endif
 183 #endif
 184
 185 #if defined CALC_ENERGIES || defined LJ_POT_SWITCH
 186                 /* Masking should be done after force switching,
 187                  * but before potential switching.
 188                  */
 189                 /* Need to zero the interaction if there should be exclusion. */
 190                 VLJ = VLJ * interact;
 191 #endif
 192
 193 #ifdef LJ_POT_SWITCH
 194                 {
 195                     real sw, dsw;
 196
 197                     sw  = 1.0 + (swV3 + (swV4 + swV5 * rsw) * rsw) * rsw * rsw * rsw;
 198                     dsw = (swF2 + (swF3 + swF4 * rsw) * rsw) * rsw * rsw;
 199
 200                     frLJ = frLJ * sw - r * VLJ * dsw;
 201                     VLJ *= sw;
 202                 }
 203 #endif
 204
 205 #ifdef LJ_EWALD
 206                 {
 207                     real c6grid, rinvsix_nm, cr2, expmcr2, poly;
 208 #    ifdef CALC_ENERGIES
 209                     real sh_mask;
 210 #    endif
 211
 212 #    ifdef LJ_EWALD_COMB_GEOM
 213                     c6grid = ljc[type[ai] * 2] * ljc[type[aj] * 2];
 214 #    elif defined LJ_EWALD_COMB_LB
 215                     {
 216                         real sigma, sigma2, epsilon;
 217
 218                         /* These sigma and epsilon are scaled to give 6*C6 */
 219                         sigma   = ljc[type[ai] * 2] + ljc[type[aj] * 2];
 220                         epsilon = ljc[type[ai] * 2 + 1] * ljc[type[aj] * 2 + 1];
 221
 222                         sigma2 = sigma * sigma;
 223                         c6grid = epsilon * sigma2 * sigma2 * sigma2;
 224                     }
 225 #    else
 226 #        error "No LJ Ewald combination rule defined"
 227 #    endif
 228
 229 #    ifdef CHECK_EXCLS
 230                     /* Recalculate rinvsix without exclusion mask */
 231                     rinvsix_nm = rinvsq * rinvsq * rinvsq;
 232 #    else
 233                     rinvsix_nm = rinvsix;
 234 #    endif
 235                     cr2 = lje_coeff2 * rsq;
 236 #    if GMX_DOUBLE
 237                     expmcr2 = exp(-cr2);
 238 #    else
 239                     expmcr2    = expf(-cr2);
 240 #    endif
 241                     poly = 1 + cr2 + 0.5 * cr2 * cr2;
 242
 243                     /* Subtract the grid force from the total LJ force */
 244                     frLJ += c6grid * (rinvsix_nm - expmcr2 * (rinvsix_nm * poly + lje_coeff6_6));
 245 #    ifdef CALC_ENERGIES
 246                     /* Shift should only be applied to real LJ pairs */
 247                     sh_mask = lje_vc * interact;
 248
 249                     VLJ += c6grid / 6 * (rinvsix_nm * (1 - expmcr2 * poly) + sh_mask);
 250 #    endif
 251                 }
 252 #endif /* LJ_EWALD */
 253
 254 #ifdef VDW_CUTOFF_CHECK
 255                 /* Mask for VdW cut-off shorter than Coulomb cut-off */
 256                 {
 257                     real skipmask_rvdw;
 258
 259                     skipmask_rvdw = (rsq < rvdw2) ? 1.0 : 0.0;
 260                     frLJ *= skipmask_rvdw;
 261 #    ifdef CALC_ENERGIES
 262                     VLJ *= skipmask_rvdw;
 263 #    endif
 264                 }
 265 #else
 266 #    if defined CALC_ENERGIES
 267                 /* Need to zero the interaction if r >= rcut */
 268                 VLJ = VLJ * skipmask;
 269                 /* 1 more flop for LJ energy */
 270 #    endif
 271 #endif /* VDW_CUTOFF_CHECK */
 272
 273
 274 #ifdef CALC_ENERGIES
 275 #    ifdef ENERGY_GROUPS
 276                 Vvdw[egp_sh_i[i] + ((egp_cj >> (nbatParams.neg_2log * j)) & egp_mask)] += VLJ;
 277 #    else
 278                 Vvdw_ci += VLJ;
 279                 /* 1 flop for LJ energy addition */
 280 #    endif
 281 #endif
 282             }
 283
 284 #ifdef CALC_COULOMB
 285             /* Enforce the cut-off and perhaps exclusions. In
 286              * those cases, rinv is zero because of skipmask,
 287              * but fcoul and vcoul will later be non-zero (in
 288              * both RF and table cases) because of the
 289              * contributions that do not depend on rinv. These
 290              * contributions cannot be allowed to accumulate
 291              * to the force and potential, and the easiest way
 292              * to do this is to zero the charges in
 293              * advance. */
 294             qq = skipmask * qi[i] * q[aj];
 295
 296 #    ifdef CALC_COUL_RF
 297             fcoul = qq * (interact * rinv * rinvsq - k_rf2);
 298             /* 4 flops for RF force */
 299 #        ifdef CALC_ENERGIES
 300             vcoul = qq * (interact * rinv + k_rf * rsq - c_rf);
 301             /* 4 flops for RF energy */
 302 #        endif
 303 #    endif
 304
 305 #    ifdef CALC_COUL_TAB
 306             rs   = rsq * rinv * tab_coul_scale;
 307             ri   = int(rs);
 308             frac = rs - static_cast<real>(ri);
 309 #        if !GMX_DOUBLE
 310             /* fexcl = F_i + frac * (F_(i+1)-F_i) */
 311             fexcl = tab_coul_FDV0[ri * 4] + frac * tab_coul_FDV0[ri * 4 + 1];
 312 #        else
 313             /* fexcl = (1-frac) * F_i + frac * F_(i+1) */
 314             fexcl = (1 - frac) * tab_coul_F[ri] + frac * tab_coul_F[ri + 1];
 315 #        endif
 316             fcoul = interact * rinvsq - fexcl;
 317             /* 7 flops for float 1/r-table force */
 318 #        ifdef CALC_ENERGIES
 319 #            if !GMX_DOUBLE
 320             vcoul = qq
 321                     * (interact * (rinv - ic->sh_ewald)
 322                        - (tab_coul_FDV0[ri * 4 + 2] - halfsp * frac * (tab_coul_FDV0[ri * 4] + fexcl)));
 323             /* 7 flops for float 1/r-table energy (8 with excls) */
 324 #            else
 325             vcoul = qq
 326                     * (interact * (rinv - ic->sh_ewald)
 327                        - (tab_coul_V[ri] - halfsp * frac * (tab_coul_F[ri] + fexcl)));
 328 #            endif
 329 #        endif
 330             fcoul *= qq * rinv;
 331 #    endif
 332
 333 #    ifdef CALC_ENERGIES
 334 #        ifdef ENERGY_GROUPS
 335             Vc[egp_sh_i[i] + ((egp_cj >> (nbatParams.neg_2log * j)) & egp_mask)] += vcoul;
 336 #        else
 337             Vc_ci += vcoul;
 338             /* 1 flop for Coulomb energy addition */
 339 #        endif
 340 #    endif
 341 #endif
 342
 343 #ifdef CALC_COULOMB
 344 #    ifdef HALF_LJ
 345             if (i < UNROLLI / 2)
 346 #    endif
 347             {
 348                 fscal = frLJ * rinvsq + fcoul;
 349                 /* 2 flops for scalar LJ+Coulomb force */
 350             }
 351 #    ifdef HALF_LJ
 352             else
 353             {
 354                 fscal = fcoul;
 355             }
 356 #    endif
 357 #else
 358             fscal = frLJ * rinvsq;
 359 #endif
 360             fx = fscal * dx;
 361             fy = fscal * dy;
 362             fz = fscal * dz;
 363
 364             /* Increment i-atom force */
 365             fi[i * FI_STRIDE + XX] += fx;
 366             fi[i * FI_STRIDE + YY] += fy;
 367             fi[i * FI_STRIDE + ZZ] += fz;
 368             /* Decrement j-atom force */
 369             f[aj * F_STRIDE + XX] -= fx;
 370             f[aj * F_STRIDE + YY] -= fy;
 371             f[aj * F_STRIDE + ZZ] -= fz;
 372             /* 9 flops for force addition */
 373         }
 374     }
 375 }
 376
 377 #undef interact
 378 #undef EXCL_FORCES