src/gromacs/ewald/pme-gather.cpp

   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-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2017, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source 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 #include "gmxpre.h"
  39
  40 #include "pme-gather.h"
  41
  42 #include "gromacs/math/vec.h"
  43 #include "gromacs/simd/simd.h"
  44 #include "gromacs/utility/basedefinitions.h"
  45 #include "gromacs/utility/smalloc.h"
  46
  47 #include "pme-internal.h"
  48 #include "pme-simd.h"
  49 #include "pme-spline-work.h"
  50
  51 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
  52
  53 #define DO_FSPLINE(order)                      \
  54     for (int ithx = 0; (ithx < order); ithx++)              \
  55     {                                              \
  56         const int  index_x = (i0+ithx)*pny*pnz;               \
  57         const real tx      = thx[ithx];                             \
  58         const real dx      = dthx[ithx];                                      \
  59                                                \
  60         for (int ithy = 0; (ithy < order); ithy++)          \
  61         {                                          \
  62             const int  index_xy = index_x+(j0+ithy)*pnz;      \
  63             const real ty       = thy[ithy];                            \
  64             const real dy       = dthy[ithy];                                     \
  65             real       fxy1     = 0, fz1 = 0;                                 \
  66                                                \
  67             for (int ithz = 0; (ithz < order); ithz++)      \
  68             {                                      \
  69                 const real gval  = grid[index_xy+(k0+ithz)];  \
  70                 fxy1 += thz[ithz]*gval;            \
  71                 fz1  += dthz[ithz]*gval;           \
  72             }                                      \
  73             fx += dx*ty*fxy1;                      \
  74             fy += tx*dy*fxy1;                      \
  75             fz += tx*ty*fz1;                       \
  76         }                                          \
  77     }
  78
  79
  80 void gather_f_bsplines(struct gmx_pme_t *pme, real *grid,
  81                        gmx_bool bClearF, pme_atomcomm_t *atc,
  82                        splinedata_t *spline,
  83                        real scale)
  84 {
  85     /* sum forces for local particles */
  86
  87 #ifdef PME_SIMD4_SPREAD_GATHER
  88     // cppcheck-suppress unreadVariable cppcheck seems not to analyze code from pme-simd4.h
  89     struct pme_spline_work *work = pme->spline_work;
  90 #ifndef PME_SIMD4_UNALIGNED
  91     GMX_ALIGNED(real, GMX_SIMD4_WIDTH)  thz_aligned[GMX_SIMD4_WIDTH*2];
  92     GMX_ALIGNED(real, GMX_SIMD4_WIDTH)  dthz_aligned[GMX_SIMD4_WIDTH*2];
  93 #endif
  94 #endif
  95
  96     const int  order = pme->pme_order;
  97     const int  nx    = pme->nkx;
  98     const int  ny    = pme->nky;
  99     const int  nz    = pme->nkz;
 100     const int  pny   = pme->pmegrid_ny;
 101     const int  pnz   = pme->pmegrid_nz;
 102
 103     const real rxx   = pme->recipbox[XX][XX];
 104     const real ryx   = pme->recipbox[YY][XX];
 105     const real ryy   = pme->recipbox[YY][YY];
 106     const real rzx   = pme->recipbox[ZZ][XX];
 107     const real rzy   = pme->recipbox[ZZ][YY];
 108     const real rzz   = pme->recipbox[ZZ][ZZ];
 109
 110     for (int nn = 0; nn < spline->n; nn++)
 111     {
 112         const int  n           = spline->ind[nn];
 113         const real coefficient = scale*atc->coefficient[n];
 114
 115         if (bClearF)
 116         {
 117             atc->f[n][XX] = 0;
 118             atc->f[n][YY] = 0;
 119             atc->f[n][ZZ] = 0;
 120         }
 121         if (coefficient != 0)
 122         {
 123             real       fx     = 0;
 124             real       fy     = 0;
 125             real       fz     = 0;
 126             const int* idxptr = atc->idx[n];
 127             const int  norder = nn*order;
 128
 129             const int  i0   = idxptr[XX];
 130             const int  j0   = idxptr[YY];
 131             const int  k0   = idxptr[ZZ];
 132
 133             /* Pointer arithmetic alert, next six statements */
 134             const real* thx  = spline->theta[XX] + norder;
 135             const real* thy  = spline->theta[YY] + norder;
 136             const real* thz  = spline->theta[ZZ] + norder;
 137             const real* dthx = spline->dtheta[XX] + norder;
 138             const real* dthy = spline->dtheta[YY] + norder;
 139             const real* dthz = spline->dtheta[ZZ] + norder;
 140
 141             switch (order)
 142             {
 143                 case 4:
 144 #ifdef PME_SIMD4_SPREAD_GATHER
 145 #ifdef PME_SIMD4_UNALIGNED
 146 #define PME_GATHER_F_SIMD4_ORDER4
 147 #else
 148 #define PME_GATHER_F_SIMD4_ALIGNED
 149 #define PME_ORDER 4
 150 #endif
 151 #include "pme-simd4.h"
 152 #else
 153                     DO_FSPLINE(4);
 154 #endif
 155                     break;
 156                 case 5:
 157 #ifdef PME_SIMD4_SPREAD_GATHER
 158 #define PME_GATHER_F_SIMD4_ALIGNED
 159 #define PME_ORDER 5
 160 #include "pme-simd4.h"
 161 #else
 162                     DO_FSPLINE(5);
 163 #endif
 164                     break;
 165                 default:
 166                     DO_FSPLINE(order);
 167                     break;
 168             }
 169
 170             atc->f[n][XX] += -coefficient*( fx*nx*rxx );
 171             atc->f[n][YY] += -coefficient*( fx*nx*ryx + fy*ny*ryy );
 172             atc->f[n][ZZ] += -coefficient*( fx*nx*rzx + fy*ny*rzy + fz*nz*rzz );
 173         }
 174     }
 175     /* Since the energy and not forces are interpolated
 176      * the net force might not be exactly zero.
 177      * This can be solved by also interpolating F, but
 178      * that comes at a cost.
 179      * A better hack is to remove the net force every
 180      * step, but that must be done at a higher level
 181      * since this routine doesn't see all atoms if running
 182      * in parallel. Don't know how important it is?  EL 990726
 183      */
 184 }
 185
 186
 187 real gather_energy_bsplines(struct gmx_pme_t *pme, real *grid,
 188                             pme_atomcomm_t *atc)
 189 {
 190     splinedata_t *spline;
 191     int     n, ithx, ithy, ithz, i0, j0, k0;
 192     int     index_x, index_xy;
 193     int *   idxptr;
 194     real    energy, pot, tx, ty, coefficient, gval;
 195     real    *thx, *thy, *thz;
 196     int     norder;
 197     int     order;
 198
 199     spline = &atc->spline[0];
 200
 201     order = pme->pme_order;
 202
 203     energy = 0;
 204     for (n = 0; (n < atc->n); n++)
 205     {
 206         coefficient      = atc->coefficient[n];
 207
 208         if (coefficient != 0)
 209         {
 210             idxptr = atc->idx[n];
 211             norder = n*order;
 212
 213             i0   = idxptr[XX];
 214             j0   = idxptr[YY];
 215             k0   = idxptr[ZZ];
 216
 217             /* Pointer arithmetic alert, next three statements */
 218             thx  = spline->theta[XX] + norder;
 219             thy  = spline->theta[YY] + norder;
 220             thz  = spline->theta[ZZ] + norder;
 221
 222             pot = 0;
 223             for (ithx = 0; (ithx < order); ithx++)
 224             {
 225                 index_x = (i0+ithx)*pme->pmegrid_ny*pme->pmegrid_nz;
 226                 tx      = thx[ithx];
 227
 228                 for (ithy = 0; (ithy < order); ithy++)
 229                 {
 230                     index_xy = index_x+(j0+ithy)*pme->pmegrid_nz;
 231                     ty       = thy[ithy];
 232
 233                     for (ithz = 0; (ithz < order); ithz++)
 234                     {
 235                         gval  = grid[index_xy+(k0+ithz)];
 236                         pot  += tx*ty*thz[ithz]*gval;
 237                     }
 238
 239                 }
 240             }
 241
 242             energy += pot*coefficient;
 243         }
 244     }
 245
 246     return energy;
 247 }