src/mdlib/nbnxn_kernels/nbnxn_kernel_ref_inner.h

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