src/mdlib/nbnxn_kernels/nbnxn_kernel_ref_outer.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,2013, 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 #define UNROLLI    NBNXN_CPU_CLUSTER_I_SIZE
  39 #define UNROLLJ    NBNXN_CPU_CLUSTER_I_SIZE
  40
  41 /* We could use nbat->xstride and nbat->fstride, but macros might be faster */
  42 #define X_STRIDE   3
  43 #define F_STRIDE   3
  44 /* Local i-atom buffer strides */
  45 #define XI_STRIDE  3
  46 #define FI_STRIDE  3
  47
  48
  49 /* All functionality defines are set here, except for:
  50  * CALC_ENERGIES, ENERGY_GROUPS which are defined before.
  51  * CHECK_EXCLS, which is set just before including the inner loop contents.
  52  */
  53
  54 /* We always calculate shift forces, because it's cheap anyhow */
  55 #define CALC_SHIFTFORCES
  56
  57 #ifdef CALC_COUL_RF
  58 #define NBK_FUNC_NAME(x,y) x##_rf_##y
  59 #endif
  60 #ifdef CALC_COUL_TAB
  61 #ifndef VDW_CUTOFF_CHECK
  62 #define NBK_FUNC_NAME(x,y) x##_tab_##y
  63 #else
  64 #define NBK_FUNC_NAME(x,y) x##_tab_twin_##y
  65 #endif
  66 #endif
  67
  68 static void
  69 #ifndef CALC_ENERGIES
  70 NBK_FUNC_NAME(nbnxn_kernel_ref,noener)
  71 #else
  72 #ifndef ENERGY_GROUPS
  73 NBK_FUNC_NAME(nbnxn_kernel_ref,ener)
  74 #else
  75 NBK_FUNC_NAME(nbnxn_kernel_ref,energrp)
  76 #endif
  77 #endif
  78 #undef NBK_FUNC_NAME
  79                             (const nbnxn_pairlist_t     *nbl,
  80                              const nbnxn_atomdata_t     *nbat,
  81                              const interaction_const_t  *ic,
  82                              rvec                       *shift_vec,
  83                              real                       *f
  84 #ifdef CALC_SHIFTFORCES
  85                              ,
  86                              real                       *fshift
  87 #endif
  88 #ifdef CALC_ENERGIES
  89                              ,
  90                              real                       *Vvdw,
  91                              real                       *Vc
  92 #endif
  93                             )
  94 {
  95     const nbnxn_ci_t   *nbln;
  96     const nbnxn_cj_t   *l_cj;
  97     const int          *type;
  98     const real         *q;
  99     const real         *shiftvec;
 100     const real         *x;
 101     const real         *nbfp;
 102     real       rcut2;
 103 #ifdef VDW_CUTOFF_CHECK
 104     real       rvdw2;
 105 #endif
 106     int        ntype2;
 107     real       facel;
 108     real       *nbfp_i;
 109     int        n,ci,ci_sh;
 110     int        ish,ishf;
 111     gmx_bool   do_LJ,half_LJ,do_coul;
 112     int        cjind0,cjind1,cjind;
 113     int        ip,jp;
 114
 115     real       xi[UNROLLI*XI_STRIDE];
 116     real       fi[UNROLLI*FI_STRIDE];
 117     real       qi[UNROLLI];
 118
 119 #ifdef CALC_ENERGIES
 120 #ifndef ENERGY_GROUPS
 121
 122     real       Vvdw_ci,Vc_ci;
 123 #else
 124     int        egp_mask;
 125     int        egp_sh_i[UNROLLI];
 126 #endif
 127     real       sh_invrc6;
 128 #endif
 129
 130 #ifdef CALC_COUL_RF
 131     real       k_rf2;
 132 #ifdef CALC_ENERGIES
 133     real       k_rf,c_rf;
 134 #endif
 135 #endif
 136 #ifdef CALC_COUL_TAB
 137     real       tabscale;
 138 #ifdef CALC_ENERGIES
 139     real       halfsp;
 140 #endif
 141 #ifndef GMX_DOUBLE
 142     const real *tab_coul_FDV0;
 143 #else
 144     const real *tab_coul_F;
 145     const real *tab_coul_V;
 146 #endif
 147 #endif
 148
 149     int ninner;
 150
 151 #ifdef COUNT_PAIRS
 152     int npair=0;
 153 #endif
 154
 155 #ifdef CALC_ENERGIES
 156     sh_invrc6 = ic->sh_invrc6;
 157 #endif
 158
 159 #ifdef CALC_COUL_RF
 160     k_rf2 = 2*ic->k_rf;
 161 #ifdef CALC_ENERGIES
 162     k_rf = ic->k_rf;
 163     c_rf = ic->c_rf;
 164 #endif
 165 #endif
 166 #ifdef CALC_COUL_TAB
 167     tabscale = ic->tabq_scale;
 168 #ifdef CALC_ENERGIES
 169     halfsp = 0.5/ic->tabq_scale;
 170 #endif
 171
 172 #ifndef GMX_DOUBLE
 173     tab_coul_FDV0 = ic->tabq_coul_FDV0;
 174 #else
 175     tab_coul_F    = ic->tabq_coul_F;
 176     tab_coul_V    = ic->tabq_coul_V;
 177 #endif
 178 #endif
 179
 180 #ifdef ENERGY_GROUPS
 181     egp_mask = (1<<nbat->neg_2log) - 1;
 182 #endif
 183
 184
 185     rcut2               = ic->rcoulomb*ic->rcoulomb;
 186 #ifdef VDW_CUTOFF_CHECK
 187     rvdw2               = ic->rvdw*ic->rvdw;
 188 #endif
 189
 190     ntype2              = nbat->ntype*2;
 191     nbfp                = nbat->nbfp;
 192     q                   = nbat->q;
 193     type                = nbat->type;
 194     facel               = ic->epsfac;
 195     shiftvec            = shift_vec[0];
 196     x                   = nbat->x;
 197
 198     l_cj = nbl->cj;
 199
 200     ninner = 0;
 201     for(n=0; n<nbl->nci; n++)
 202     {
 203         int i,d;
 204
 205         nbln = &nbl->ci[n];
 206
 207         ish              = (nbln->shift & NBNXN_CI_SHIFT);
 208         /* x, f and fshift are assumed to be stored with stride 3 */
 209         ishf             = ish*DIM;
 210         cjind0           = nbln->cj_ind_start;
 211         cjind1           = nbln->cj_ind_end;
 212         /* Currently only works super-cells equal to sub-cells */
 213         ci               = nbln->ci;
 214         ci_sh            = (ish == CENTRAL ? ci : -1);
 215
 216         /* We have 5 LJ/C combinations, but use only three inner loops,
 217          * as the other combinations are unlikely and/or not much faster:
 218          * inner half-LJ + C for half-LJ + C / no-LJ + C
 219          * inner LJ + C      for full-LJ + C
 220          * inner LJ          for full-LJ + no-C / half-LJ + no-C
 221          */
 222         do_LJ   = (nbln->shift & NBNXN_CI_DO_LJ(0));
 223         do_coul = (nbln->shift & NBNXN_CI_DO_COUL(0));
 224         half_LJ = ((nbln->shift & NBNXN_CI_HALF_LJ(0)) || !do_LJ) && do_coul;
 225
 226 #ifdef CALC_ENERGIES
 227 #ifndef ENERGY_GROUPS
 228         Vvdw_ci = 0;
 229         Vc_ci   = 0;
 230 #else
 231         for(i=0; i<UNROLLI; i++)
 232         {
 233             egp_sh_i[i] = ((nbat->energrp[ci]>>(i*nbat->neg_2log)) & egp_mask)*nbat->nenergrp;
 234         }
 235 #endif
 236 #endif
 237
 238         for(i=0; i<UNROLLI; i++)
 239         {
 240             for(d=0; d<DIM; d++)
 241             {
 242                 xi[i*XI_STRIDE+d] = x[(ci*UNROLLI+i)*X_STRIDE+d] + shiftvec[ishf+d];
 243                 fi[i*FI_STRIDE+d] = 0;
 244             }
 245         }
 246
 247         if (do_coul)
 248         {
 249 #ifdef CALC_ENERGIES
 250             real Vc_sub_self;
 251
 252 #ifdef CALC_COUL_RF
 253             Vc_sub_self = 0.5*c_rf;
 254 #endif
 255 #ifdef CALC_COUL_TAB
 256 #ifdef GMX_DOUBLE
 257             Vc_sub_self = 0.5*tab_coul_V[0];
 258 #else
 259             Vc_sub_self = 0.5*tab_coul_FDV0[2];
 260 #endif
 261 #endif
 262 #endif
 263
 264             for(i=0; i<UNROLLI; i++)
 265             {
 266                 qi[i] = facel*q[ci*UNROLLI+i];
 267
 268 #ifdef CALC_ENERGIES
 269                 if (l_cj[nbln->cj_ind_start].cj == ci_sh)
 270                 {
 271 #ifdef ENERGY_GROUPS
 272                     Vc[egp_sh_i[i]+((nbat->energrp[ci]>>(i*nbat->neg_2log)) & egp_mask)]
 273 #else
 274                     Vc[0]
 275 #endif
 276                         -= qi[i]*q[ci*UNROLLI+i]*Vc_sub_self;
 277                 }
 278 #endif
 279             }
 280         }
 281
 282         cjind = cjind0;
 283         while (cjind < cjind1 && nbl->cj[cjind].excl != 0xffff)
 284         {
 285 #define CHECK_EXCLS
 286             if (half_LJ)
 287             {
 288 #define CALC_COULOMB
 289 #define HALF_LJ
 290 #include "nbnxn_kernel_ref_inner.h"
 291 #undef HALF_LJ
 292 #undef CALC_COULOMB
 293             }
 294             /* cppcheck-suppress duplicateBranch */
 295             else if (do_coul)
 296             {
 297 #define CALC_COULOMB
 298 #include "nbnxn_kernel_ref_inner.h"
 299 #undef CALC_COULOMB
 300             }
 301             else
 302             {
 303 #include "nbnxn_kernel_ref_inner.h"
 304             }
 305 #undef CHECK_EXCLS
 306             cjind++;
 307         }
 308
 309         for(; (cjind<cjind1); cjind++)
 310         {
 311             if (half_LJ)
 312             {
 313 #define CALC_COULOMB
 314 #define HALF_LJ
 315 #include "nbnxn_kernel_ref_inner.h"
 316 #undef HALF_LJ
 317 #undef CALC_COULOMB
 318             }
 319             /* cppcheck-suppress duplicateBranch */
 320             else if (do_coul)
 321             {
 322 #define CALC_COULOMB
 323 #include "nbnxn_kernel_ref_inner.h"
 324 #undef CALC_COULOMB
 325             }
 326             else
 327             {
 328 #include "nbnxn_kernel_ref_inner.h"
 329             }
 330         }
 331         ninner += cjind1 - cjind0;
 332
 333         /* Add accumulated i-forces to the force array */
 334         for(i=0; i<UNROLLI; i++)
 335         {
 336             for(d=0; d<DIM; d++)
 337             {
 338                 f[(ci*UNROLLI+i)*F_STRIDE+d] += fi[i*FI_STRIDE+d];
 339             }
 340         }
 341 #ifdef CALC_SHIFTFORCES
 342         if (fshift != NULL)
 343         {
 344             /* Add i forces to shifted force list */
 345             for(i=0; i<UNROLLI; i++)
 346             {
 347                 for(d=0; d<DIM; d++)
 348                 {
 349                     fshift[ishf+d] += fi[i*FI_STRIDE+d];
 350                 }
 351             }
 352         }
 353 #endif
 354
 355 #ifdef CALC_ENERGIES
 356 #ifndef ENERGY_GROUPS
 357         *Vvdw += Vvdw_ci;
 358         *Vc   += Vc_ci;
 359 #endif
 360 #endif
 361         }
 362
 363 #ifdef COUNT_PAIRS
 364     printf("atom pairs %d\n",npair);
 365 #endif
 366 }
 367
 368 #undef CALC_SHIFTFORCES
 369
 370 #undef X_STRIDE
 371 #undef F_STRIDE
 372 #undef XI_STRIDE
 373 #undef FI_STRIDE
 374
 375 #undef UNROLLI
 376 #undef UNROLLJ