src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref_outer.h

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   9  * Copyright (c) 2001-2009, The GROMACS Development Team
  10  *
  11  * Gromacs is a library for molecular simulation and trajectory analysis,
  12  * written by Erik Lindahl, David van der Spoel, Berk Hess, and others - for
  13  * a full list of developers and information, check out http://www.gromacs.org
  14  *
  15  * This program is free software; you can redistribute it and/or modify it under
  16  * the terms of the GNU Lesser General Public License as published by the Free
  17  * Software Foundation; either version 2 of the License, or (at your option) any
  18  * later version.
  19  * As a special exception, you may use this file as part of a free software
  20  * library without restriction.  Specifically, if other files instantiate
  21  * templates or use macros or inline functions from this file, or you compile
  22  * this file and link it with other files to produce an executable, this
  23  * file does not by itself cause the resulting executable to be covered by
  24  * the GNU Lesser General Public License.
  25  *
  26  * In plain-speak: do not worry about classes/macros/templates either - only
  27  * changes to the library have to be LGPL, not an application linking with it.
  28  *
  29  * To help fund GROMACS development, we humbly ask that you cite
  30  * the papers people have written on it - you can find them on the website!
  31  */
  32
  33 #define UNROLLI    NBNXN_CPU_CLUSTER_I_SIZE
  34 #define UNROLLJ    NBNXN_CPU_CLUSTER_I_SIZE
  35
  36 /* We could use nbat->xstride and nbat->fstride, but macros might be faster */
  37 #define X_STRIDE   3
  38 #define F_STRIDE   3
  39 /* Local i-atom buffer strides */
  40 #define XI_STRIDE  3
  41 #define FI_STRIDE  3
  42
  43
  44 /* All functionality defines are set here, except for:
  45  * CALC_ENERGIES, ENERGY_GROUPS which are defined before.
  46  * CHECK_EXCLS, which is set just before including the inner loop contents.
  47  */
  48
  49 /* We always calculate shift forces, because it's cheap anyhow */
  50 #define CALC_SHIFTFORCES
  51
  52 #ifdef CALC_COUL_RF
  53 #define NBK_FUNC_NAME(base, ene) base ## _rf_ ## ene
  54 #endif
  55 #ifdef CALC_COUL_TAB
  56 #ifndef VDW_CUTOFF_CHECK
  57 #define NBK_FUNC_NAME(base, ene) base ## _tab_ ## ene
  58 #else
  59 #define NBK_FUNC_NAME(base, ene) base ## _tab_twin_ ## ene
  60 #endif
  61 #endif
  62
  63 static void
  64 #ifndef CALC_ENERGIES
  65 NBK_FUNC_NAME(nbnxn_kernel_ref, noener)
  66 #else
  67 #ifndef ENERGY_GROUPS
  68 NBK_FUNC_NAME(nbnxn_kernel_ref, ener)
  69 #else
  70 NBK_FUNC_NAME(nbnxn_kernel_ref, energrp)
  71 #endif
  72 #endif
  73 #undef NBK_FUNC_NAME
  74 (const nbnxn_pairlist_t     *nbl,
  75  const nbnxn_atomdata_t     *nbat,
  76  const interaction_const_t  *ic,
  77  rvec                       *shift_vec,
  78  real                       *f
  79 #ifdef CALC_SHIFTFORCES
  80  ,
  81  real                       *fshift
  82 #endif
  83 #ifdef CALC_ENERGIES
  84  ,
  85  real                       *Vvdw,
  86  real                       *Vc
  87 #endif
  88 )
  89 {
  90     const nbnxn_ci_t   *nbln;
  91     const nbnxn_cj_t   *l_cj;
  92     const int          *type;
  93     const real         *q;
  94     const real         *shiftvec;
  95     const real         *x;
  96     const real         *nbfp;
  97     real                rcut2;
  98 #ifdef VDW_CUTOFF_CHECK
  99     real                rvdw2;
 100 #endif
 101     int                 ntype2;
 102     real                facel;
 103     real               *nbfp_i;
 104     int                 n, ci, ci_sh;
 105     int                 ish, ishf;
 106     gmx_bool            do_LJ, half_LJ, do_coul;
 107     int                 cjind0, cjind1, cjind;
 108     int                 ip, jp;
 109
 110     real                xi[UNROLLI*XI_STRIDE];
 111     real                fi[UNROLLI*FI_STRIDE];
 112     real                qi[UNROLLI];
 113
 114 #ifdef CALC_ENERGIES
 115 #ifndef ENERGY_GROUPS
 116
 117     real       Vvdw_ci, Vc_ci;
 118 #else
 119     int        egp_mask;
 120     int        egp_sh_i[UNROLLI];
 121 #endif
 122     real       sh_invrc6;
 123 #endif
 124
 125 #ifdef CALC_COUL_RF
 126     real       k_rf2;
 127 #ifdef CALC_ENERGIES
 128     real       k_rf, c_rf;
 129 #endif
 130 #endif
 131 #ifdef CALC_COUL_TAB
 132     real       tabscale;
 133 #ifdef CALC_ENERGIES
 134     real       halfsp;
 135 #endif
 136 #ifndef GMX_DOUBLE
 137     const real *tab_coul_FDV0;
 138 #else
 139     const real *tab_coul_F;
 140     const real *tab_coul_V;
 141 #endif
 142 #endif
 143
 144     int ninner;
 145
 146 #ifdef COUNT_PAIRS
 147     int npair = 0;
 148 #endif
 149
 150 #ifdef CALC_ENERGIES
 151     sh_invrc6 = ic->sh_invrc6;
 152 #endif
 153
 154 #ifdef CALC_COUL_RF
 155     k_rf2 = 2*ic->k_rf;
 156 #ifdef CALC_ENERGIES
 157     k_rf = ic->k_rf;
 158     c_rf = ic->c_rf;
 159 #endif
 160 #endif
 161 #ifdef CALC_COUL_TAB
 162     tabscale = ic->tabq_scale;
 163 #ifdef CALC_ENERGIES
 164     halfsp = 0.5/ic->tabq_scale;
 165 #endif
 166
 167 #ifndef GMX_DOUBLE
 168     tab_coul_FDV0 = ic->tabq_coul_FDV0;
 169 #else
 170     tab_coul_F    = ic->tabq_coul_F;
 171     tab_coul_V    = ic->tabq_coul_V;
 172 #endif
 173 #endif
 174
 175 #ifdef ENERGY_GROUPS
 176     egp_mask = (1<<nbat->neg_2log) - 1;
 177 #endif
 178
 179
 180     rcut2               = ic->rcoulomb*ic->rcoulomb;
 181 #ifdef VDW_CUTOFF_CHECK
 182     rvdw2               = ic->rvdw*ic->rvdw;
 183 #endif
 184
 185     ntype2              = nbat->ntype*2;
 186     nbfp                = nbat->nbfp;
 187     q                   = nbat->q;
 188     type                = nbat->type;
 189     facel               = ic->epsfac;
 190     shiftvec            = shift_vec[0];
 191     x                   = nbat->x;
 192
 193     l_cj = nbl->cj;
 194
 195     ninner = 0;
 196     for (n = 0; n < nbl->nci; n++)
 197     {
 198         int i, d;
 199
 200         nbln = &nbl->ci[n];
 201
 202         ish              = (nbln->shift & NBNXN_CI_SHIFT);
 203         /* x, f and fshift are assumed to be stored with stride 3 */
 204         ishf             = ish*DIM;
 205         cjind0           = nbln->cj_ind_start;
 206         cjind1           = nbln->cj_ind_end;
 207         /* Currently only works super-cells equal to sub-cells */
 208         ci               = nbln->ci;
 209         ci_sh            = (ish == CENTRAL ? ci : -1);
 210
 211         /* We have 5 LJ/C combinations, but use only three inner loops,
 212          * as the other combinations are unlikely and/or not much faster:
 213          * inner half-LJ + C for half-LJ + C / no-LJ + C
 214          * inner LJ + C      for full-LJ + C
 215          * inner LJ          for full-LJ + no-C / half-LJ + no-C
 216          */
 217         do_LJ   = (nbln->shift & NBNXN_CI_DO_LJ(0));
 218         do_coul = (nbln->shift & NBNXN_CI_DO_COUL(0));
 219         half_LJ = ((nbln->shift & NBNXN_CI_HALF_LJ(0)) || !do_LJ) && do_coul;
 220
 221 #ifdef CALC_ENERGIES
 222 #ifndef ENERGY_GROUPS
 223         Vvdw_ci = 0;
 224         Vc_ci   = 0;
 225 #else
 226         for (i = 0; i < UNROLLI; i++)
 227         {
 228             egp_sh_i[i] = ((nbat->energrp[ci]>>(i*nbat->neg_2log)) & egp_mask)*nbat->nenergrp;
 229         }
 230 #endif
 231 #endif
 232
 233         for (i = 0; i < UNROLLI; i++)
 234         {
 235             for (d = 0; d < DIM; d++)
 236             {
 237                 xi[i*XI_STRIDE+d] = x[(ci*UNROLLI+i)*X_STRIDE+d] + shiftvec[ishf+d];
 238                 fi[i*FI_STRIDE+d] = 0;
 239             }
 240         }
 241
 242         if (do_coul)
 243         {
 244 #ifdef CALC_ENERGIES
 245             real Vc_sub_self;
 246
 247 #ifdef CALC_COUL_RF
 248             Vc_sub_self = 0.5*c_rf;
 249 #endif
 250 #ifdef CALC_COUL_TAB
 251 #ifdef GMX_DOUBLE
 252             Vc_sub_self = 0.5*tab_coul_V[0];
 253 #else
 254             Vc_sub_self = 0.5*tab_coul_FDV0[2];
 255 #endif
 256 #endif
 257 #endif
 258
 259             for (i = 0; i < UNROLLI; i++)
 260             {
 261                 qi[i] = facel*q[ci*UNROLLI+i];
 262
 263 #ifdef CALC_ENERGIES
 264                 if (l_cj[nbln->cj_ind_start].cj == ci_sh)
 265                 {
 266 #ifdef ENERGY_GROUPS
 267                     Vc[egp_sh_i[i]+((nbat->energrp[ci]>>(i*nbat->neg_2log)) & egp_mask)]
 268 #else
 269                     Vc[0]
 270 #endif
 271                         -= qi[i]*q[ci*UNROLLI+i]*Vc_sub_self;
 272                 }
 273 #endif
 274             }
 275         }
 276
 277         cjind = cjind0;
 278         while (cjind < cjind1 && nbl->cj[cjind].excl != 0xffff)
 279         {
 280 #define CHECK_EXCLS
 281             if (half_LJ)
 282             {
 283 #define CALC_COULOMB
 284 #define HALF_LJ
 285 #include "nbnxn_kernel_ref_inner.h"
 286 #undef HALF_LJ
 287 #undef CALC_COULOMB
 288             }
 289             /* cppcheck-suppress duplicateBranch */
 290             else if (do_coul)
 291             {
 292 #define CALC_COULOMB
 293 #include "nbnxn_kernel_ref_inner.h"
 294 #undef CALC_COULOMB
 295             }
 296             else
 297             {
 298 #include "nbnxn_kernel_ref_inner.h"
 299             }
 300 #undef CHECK_EXCLS
 301             cjind++;
 302         }
 303
 304         for (; (cjind < cjind1); cjind++)
 305         {
 306             if (half_LJ)
 307             {
 308 #define CALC_COULOMB
 309 #define HALF_LJ
 310 #include "nbnxn_kernel_ref_inner.h"
 311 #undef HALF_LJ
 312 #undef CALC_COULOMB
 313             }
 314             /* cppcheck-suppress duplicateBranch */
 315             else if (do_coul)
 316             {
 317 #define CALC_COULOMB
 318 #include "nbnxn_kernel_ref_inner.h"
 319 #undef CALC_COULOMB
 320             }
 321             else
 322             {
 323 #include "nbnxn_kernel_ref_inner.h"
 324             }
 325         }
 326         ninner += cjind1 - cjind0;
 327
 328         /* Add accumulated i-forces to the force array */
 329         for (i = 0; i < UNROLLI; i++)
 330         {
 331             for (d = 0; d < DIM; d++)
 332             {
 333                 f[(ci*UNROLLI+i)*F_STRIDE+d] += fi[i*FI_STRIDE+d];
 334             }
 335         }
 336 #ifdef CALC_SHIFTFORCES
 337         if (fshift != NULL)
 338         {
 339             /* Add i forces to shifted force list */
 340             for (i = 0; i < UNROLLI; i++)
 341             {
 342                 for (d = 0; d < DIM; d++)
 343                 {
 344                     fshift[ishf+d] += fi[i*FI_STRIDE+d];
 345                 }
 346             }
 347         }
 348 #endif
 349
 350 #ifdef CALC_ENERGIES
 351 #ifndef ENERGY_GROUPS
 352         *Vvdw += Vvdw_ci;
 353         *Vc   += Vc_ci;
 354 #endif
 355 #endif
 356     }
 357
 358 #ifdef COUNT_PAIRS
 359     printf("atom pairs %d\n", npair);
 360 #endif
 361 }
 362
 363 #undef CALC_SHIFTFORCES
 364
 365 #undef X_STRIDE
 366 #undef F_STRIDE
 367 #undef XI_STRIDE
 368 #undef FI_STRIDE
 369
 370 #undef UNROLLI
 371 #undef UNROLLJ