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(x,y) x##_rf_##y
  54 #endif
  55 #ifdef CALC_COUL_TAB
  56 #ifndef VDW_CUTOFF_CHECK
  57 #define NBK_FUNC_NAME(x,y) x##_tab_##y
  58 #else
  59 #define NBK_FUNC_NAME(x,y) x##_tab_twin_##y
  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   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         half_LJ = (nbln->shift & NBNXN_CI_HALF_LJ(0));
 212         do_coul = (nbln->shift & NBNXN_CI_DO_COUL(0));
 213
 214 #ifdef CALC_ENERGIES
 215 #ifndef ENERGY_GROUPS
 216         Vvdw_ci = 0;
 217         Vc_ci   = 0;
 218 #else
 219         for(i=0; i<UNROLLI; i++)
 220         {
 221             egp_sh_i[i] = ((nbat->energrp[ci]>>(i*nbat->neg_2log)) & egp_mask)*nbat->nenergrp;
 222         }
 223 #endif
 224 #endif
 225
 226         for(i=0; i<UNROLLI; i++)
 227         {
 228             for(d=0; d<DIM; d++)
 229             {
 230                 xi[i*XI_STRIDE+d] = x[(ci*UNROLLI+i)*X_STRIDE+d] + shiftvec[ishf+d];
 231                 fi[i*FI_STRIDE+d] = 0;
 232             }
 233         }
 234
 235         /* With half_LJ we currently always calculate Coulomb interactions */
 236         if (do_coul || half_LJ)
 237         {
 238 #ifdef CALC_ENERGIES
 239             real Vc_sub_self;
 240
 241 #ifdef CALC_COUL_RF
 242             Vc_sub_self = 0.5*c_rf;
 243 #endif
 244 #ifdef CALC_COUL_TAB
 245 #ifdef GMX_DOUBLE
 246             Vc_sub_self = 0.5*tab_coul_V[0];
 247 #else
 248             Vc_sub_self = 0.5*tab_coul_FDV0[2];
 249 #endif
 250 #endif
 251 #endif
 252
 253             for(i=0; i<UNROLLI; i++)
 254             {
 255                 qi[i] = facel*q[ci*UNROLLI+i];
 256
 257 #ifdef CALC_ENERGIES
 258                 if (l_cj[nbln->cj_ind_start].cj == ci_sh)
 259                 {
 260 #ifdef ENERGY_GROUPS
 261                     Vc[egp_sh_i[i]+((nbat->energrp[ci]>>(i*nbat->neg_2log)) & egp_mask)]
 262 #else
 263                     Vc[0]
 264 #endif
 265                         -= qi[i]*q[ci*UNROLLI+i]*Vc_sub_self;
 266                 }
 267 #endif
 268             }
 269         }
 270
 271         cjind = cjind0;
 272         while (cjind < cjind1 && nbl->cj[cjind].excl != 0xffff)
 273         {
 274 #define CHECK_EXCLS
 275             if (half_LJ)
 276             {
 277 #define CALC_COULOMB
 278 #define HALF_LJ
 279 #include "nbnxn_kernel_ref_inner.h"
 280 #undef HALF_LJ
 281 #undef CALC_COULOMB
 282             }
 283             /* cppcheck-suppress duplicateBranch */
 284             else if (do_coul)
 285             {
 286 #define CALC_COULOMB
 287 #include "nbnxn_kernel_ref_inner.h"
 288 #undef CALC_COULOMB
 289             }
 290             else
 291             {
 292 #include "nbnxn_kernel_ref_inner.h"
 293             }
 294 #undef CHECK_EXCLS
 295             cjind++;
 296         }
 297
 298         for(; (cjind<cjind1); cjind++)
 299         {
 300             if (half_LJ)
 301             {
 302 #define CALC_COULOMB
 303 #define HALF_LJ
 304 #include "nbnxn_kernel_ref_inner.h"
 305 #undef HALF_LJ
 306 #undef CALC_COULOMB
 307             }
 308             /* cppcheck-suppress duplicateBranch */
 309             else if (do_coul)
 310             {
 311 #define CALC_COULOMB
 312 #include "nbnxn_kernel_ref_inner.h"
 313 #undef CALC_COULOMB
 314             }
 315             else
 316             {
 317 #include "nbnxn_kernel_ref_inner.h"
 318             }
 319         }
 320         ninner += cjind1 - cjind0;
 321
 322         /* Add accumulated i-forces to the force array */
 323         for(i=0; i<UNROLLI; i++)
 324         {
 325             for(d=0; d<DIM; d++)
 326             {
 327                 f[(ci*UNROLLI+i)*F_STRIDE+d] += fi[i*FI_STRIDE+d];
 328             }
 329         }
 330 #ifdef CALC_SHIFTFORCES
 331         if (fshift != NULL)
 332         {
 333             /* Add i forces to shifted force list */
 334             for(i=0; i<UNROLLI; i++)
 335             {
 336                 for(d=0; d<DIM; d++)
 337                 {
 338                     fshift[ishf+d] += fi[i*FI_STRIDE+d];
 339                 }
 340             }
 341         }
 342 #endif
 343
 344 #ifdef CALC_ENERGIES
 345 #ifndef ENERGY_GROUPS
 346         *Vvdw += Vvdw_ci;
 347         *Vc   += Vc_ci;
 348 #endif
 349 #endif
 350         }
 351
 352 #ifdef COUNT_PAIRS
 353     printf("atom pairs %d\n",npair);
 354 #endif
 355 }
 356
 357 #undef CALC_SHIFTFORCES
 358
 359 #undef X_STRIDE
 360 #undef F_STRIDE
 361 #undef XI_STRIDE
 362 #undef FI_STRIDE
 363
 364 #undef UNROLLI
 365 #undef UNROLLJ