src/gromacs/ewald/pme.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,2016,2017,2018, 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 /*! \internal \file
  38  *
  39  * \brief This file contains function definitions necessary for
  40  * computing energies and forces for the PME long-ranged part (Coulomb
  41  * and LJ).
  42  *
  43  * \author Erik Lindahl <erik@kth.se>
  44  * \author Berk Hess <hess@kth.se>
  45  * \ingroup module_ewald
  46  */
  47 /* IMPORTANT FOR DEVELOPERS:
  48  *
  49  * Triclinic pme stuff isn't entirely trivial, and we've experienced
  50  * some bugs during development (many of them due to me). To avoid
  51  * this in the future, please check the following things if you make
  52  * changes in this file:
  53  *
  54  * 1. You should obtain identical (at least to the PME precision)
  55  *    energies, forces, and virial for
  56  *    a rectangular box and a triclinic one where the z (or y) axis is
  57  *    tilted a whole box side. For instance you could use these boxes:
  58  *
  59  *    rectangular       triclinic
  60  *     2  0  0           2  0  0
  61  *     0  2  0           0  2  0
  62  *     0  0  6           2  2  6
  63  *
  64  * 2. You should check the energy conservation in a triclinic box.
  65  *
  66  * It might seem an overkill, but better safe than sorry.
  67  * /Erik 001109
  68  */
  69
  70 #include "gmxpre.h"
  71
  72 #include "pme.h"
  73
  74 #include "config.h"
  75
  76 #include <cassert>
  77 #include <cmath>
  78 #include <cstdio>
  79 #include <cstdlib>
  80 #include <cstring>
  81
  82 #include <algorithm>
  83 #include <list>
  84
  85 #include "gromacs/domdec/domdec.h"
  86 #include "gromacs/ewald/ewald-utils.h"
  87 #include "gromacs/fft/parallel_3dfft.h"
  88 #include "gromacs/fileio/pdbio.h"
  89 #include "gromacs/gmxlib/network.h"
  90 #include "gromacs/gmxlib/nrnb.h"
  91 #include "gromacs/math/gmxcomplex.h"
  92 #include "gromacs/math/invertmatrix.h"
  93 #include "gromacs/math/units.h"
  94 #include "gromacs/math/vec.h"
  95 #include "gromacs/math/vectypes.h"
  96 #include "gromacs/mdtypes/commrec.h"
  97 #include "gromacs/mdtypes/forcerec.h"
  98 #include "gromacs/mdtypes/inputrec.h"
  99 #include "gromacs/mdtypes/md_enums.h"
 100 #include "gromacs/pbcutil/pbc.h"
 101 #include "gromacs/timing/cyclecounter.h"
 102 #include "gromacs/timing/wallcycle.h"
 103 #include "gromacs/timing/walltime_accounting.h"
 104 #include "gromacs/topology/topology.h"
 105 #include "gromacs/utility/basedefinitions.h"
 106 #include "gromacs/utility/exceptions.h"
 107 #include "gromacs/utility/fatalerror.h"
 108 #include "gromacs/utility/gmxmpi.h"
 109 #include "gromacs/utility/gmxomp.h"
 110 #include "gromacs/utility/logger.h"
 111 #include "gromacs/utility/real.h"
 112 #include "gromacs/utility/smalloc.h"
 113 #include "gromacs/utility/stringutil.h"
 114 #include "gromacs/utility/unique_cptr.h"
 115
 116 #include "calculate-spline-moduli.h"
 117 #include "pme-gather.h"
 118 #include "pme-gpu-internal.h"
 119 #include "pme-grid.h"
 120 #include "pme-internal.h"
 121 #include "pme-redistribute.h"
 122 #include "pme-solve.h"
 123 #include "pme-spline-work.h"
 124 #include "pme-spread.h"
 125
 126 /*! \brief Help build a descriptive message in \c error if there are
 127  * \c errorReasons why PME on GPU is not supported.
 128  *
 129  * \returns Whether the lack of errorReasons indicate there is support. */
 130 static bool
 131 addMessageIfNotSupported(const std::list<std::string> &errorReasons,
 132                          std::string                  *error)
 133 {
 134     bool foundErrorReasons = errorReasons.empty();
 135     if (!foundErrorReasons && error)
 136     {
 137         std::string regressionTestMarker = "PME GPU does not support";
 138         // this prefix is tested for in the regression tests script gmxtest.pl
 139         *error = regressionTestMarker + ": " + gmx::joinStrings(errorReasons, "; ") + ".";
 140     }
 141     return foundErrorReasons;
 142 }
 143
 144 bool pme_gpu_supports_build(std::string *error)
 145 {
 146     std::list<std::string> errorReasons;
 147     if (GMX_DOUBLE)
 148     {
 149         errorReasons.emplace_back("double precision");
 150     }
 151     if (GMX_GPU == GMX_GPU_NONE)
 152     {
 153         errorReasons.emplace_back("non-GPU build of GROMACS");
 154     }
 155     return addMessageIfNotSupported(errorReasons, error);
 156 }
 157
 158 bool pme_gpu_supports_input(const t_inputrec &ir, const gmx_mtop_t &mtop, std::string *error)
 159 {
 160     std::list<std::string> errorReasons;
 161     if (!EEL_PME(ir.coulombtype))
 162     {
 163         errorReasons.emplace_back("systems that do not use PME for electrostatics");
 164     }
 165     if (ir.pme_order != 4)
 166     {
 167         errorReasons.emplace_back("interpolation orders other than 4");
 168     }
 169     if (ir.efep != efepNO)
 170     {
 171         if (gmx_mtop_has_perturbed_charges(mtop))
 172         {
 173             errorReasons.emplace_back("free energy calculations with perturbed charges (multiple grids)");
 174         }
 175     }
 176     if (EVDW_PME(ir.vdwtype))
 177     {
 178         errorReasons.emplace_back("Lennard-Jones PME");
 179     }
 180     if (ir.cutoff_scheme == ecutsGROUP)
 181     {
 182         errorReasons.emplace_back("group cutoff scheme");
 183     }
 184     if (!EI_DYNAMICS(ir.eI))
 185     {
 186         errorReasons.emplace_back("not a dynamical integrator");
 187     }
 188     return addMessageIfNotSupported(errorReasons, error);
 189 }
 190
 191 /*! \brief \libinternal
 192  * Finds out if PME with given inputs is possible to run on GPU.
 193  * This function is an internal final check, validating the whole PME structure on creation,
 194  * but it still duplicates the preliminary checks from the above (externally exposed) pme_gpu_supports_input() - just in case.
 195  *
 196  * \param[in]  pme          The PME structure.
 197  * \param[out] error        The error message if the input is not supported on GPU.
 198  * \returns                 True if this PME input is possible to run on GPU, false otherwise.
 199  */
 200 static bool pme_gpu_check_restrictions(const gmx_pme_t *pme, std::string *error)
 201 {
 202     std::list<std::string> errorReasons;
 203     if (pme->nnodes != 1)
 204     {
 205         errorReasons.emplace_back("PME decomposition");
 206     }
 207     if (pme->pme_order != 4)
 208     {
 209         errorReasons.emplace_back("interpolation orders other than 4");
 210     }
 211     if (pme->bFEP)
 212     {
 213         errorReasons.emplace_back("free energy calculations (multiple grids)");
 214     }
 215     if (pme->doLJ)
 216     {
 217         errorReasons.emplace_back("Lennard-Jones PME");
 218     }
 219     if (GMX_DOUBLE)
 220     {
 221         errorReasons.emplace_back("double precision");
 222     }
 223     if (GMX_GPU == GMX_GPU_NONE)
 224     {
 225         errorReasons.emplace_back("non-GPU build of GROMACS");
 226     }
 227
 228     return addMessageIfNotSupported(errorReasons, error);
 229 }
 230
 231 PmeRunMode pme_run_mode(const gmx_pme_t *pme)
 232 {
 233     GMX_ASSERT(pme != nullptr, "Expecting valid PME data pointer");
 234     return pme->runMode;
 235 }
 236
 237 gmx::PinningPolicy pme_get_pinning_policy()
 238 {
 239     return gmx::PinningPolicy::PinnedIfSupported;
 240 }
 241
 242 /*! \brief Number of bytes in a cache line.
 243  *
 244  * Must also be a multiple of the SIMD and SIMD4 register size, to
 245  * preserve alignment.
 246  */
 247 const int gmxCacheLineSize = 64;
 248
 249 //! Set up coordinate communication
 250 static void setup_coordinate_communication(pme_atomcomm_t *atc)
 251 {
 252     int nslab, n, i;
 253     int fw, bw;
 254
 255     nslab = atc->nslab;
 256
 257     n = 0;
 258     for (i = 1; i <= nslab/2; i++)
 259     {
 260         fw = (atc->nodeid + i) % nslab;
 261         bw = (atc->nodeid - i + nslab) % nslab;
 262         if (n < nslab - 1)
 263         {
 264             atc->node_dest[n] = fw;
 265             atc->node_src[n]  = bw;
 266             n++;
 267         }
 268         if (n < nslab - 1)
 269         {
 270             atc->node_dest[n] = bw;
 271             atc->node_src[n]  = fw;
 272             n++;
 273         }
 274     }
 275 }
 276
 277 /*! \brief Round \p n up to the next multiple of \p f */
 278 static int mult_up(int n, int f)
 279 {
 280     return ((n + f - 1)/f)*f;
 281 }
 282
 283 /*! \brief Return estimate of the load imbalance from the PME grid not being a good match for the number of PME ranks */
 284 static double estimate_pme_load_imbalance(struct gmx_pme_t *pme)
 285 {
 286     int    nma, nmi;
 287     double n1, n2, n3;
 288
 289     nma = pme->nnodes_major;
 290     nmi = pme->nnodes_minor;
 291
 292     n1 = mult_up(pme->nkx, nma)*mult_up(pme->nky, nmi)*pme->nkz;
 293     n2 = mult_up(pme->nkx, nma)*mult_up(pme->nkz, nmi)*pme->nky;
 294     n3 = mult_up(pme->nky, nma)*mult_up(pme->nkz, nmi)*pme->nkx;
 295
 296     /* pme_solve is roughly double the cost of an fft */
 297
 298     return (n1 + n2 + 3*n3)/static_cast<double>(6*pme->nkx*pme->nky*pme->nkz);
 299 }
 300
 301 /*! \brief Initialize atom communication data structure */
 302 static void init_atomcomm(struct gmx_pme_t *pme, pme_atomcomm_t *atc,
 303                           int dimind, gmx_bool bSpread)
 304 {
 305     int thread;
 306
 307     atc->dimind    = dimind;
 308     atc->nslab     = 1;
 309     atc->nodeid    = 0;
 310     atc->pd_nalloc = 0;
 311 #if GMX_MPI
 312     if (pme->nnodes > 1)
 313     {
 314         atc->mpi_comm = pme->mpi_comm_d[dimind];
 315         MPI_Comm_size(atc->mpi_comm, &atc->nslab);
 316         MPI_Comm_rank(atc->mpi_comm, &atc->nodeid);
 317     }
 318     if (debug)
 319     {
 320         fprintf(debug, "For PME atom communication in dimind %d: nslab %d rank %d\n", atc->dimind, atc->nslab, atc->nodeid);
 321     }
 322 #endif
 323
 324     atc->bSpread   = bSpread;
 325     atc->pme_order = pme->pme_order;
 326
 327     if (atc->nslab > 1)
 328     {
 329         snew(atc->node_dest, atc->nslab);
 330         snew(atc->node_src, atc->nslab);
 331         setup_coordinate_communication(atc);
 332
 333         snew(atc->count_thread, pme->nthread);
 334         for (thread = 0; thread < pme->nthread; thread++)
 335         {
 336             snew(atc->count_thread[thread], atc->nslab);
 337         }
 338         atc->count = atc->count_thread[0];
 339         snew(atc->rcount, atc->nslab);
 340         snew(atc->buf_index, atc->nslab);
 341     }
 342
 343     atc->nthread = pme->nthread;
 344     if (atc->nthread > 1)
 345     {
 346         snew(atc->thread_plist, atc->nthread);
 347     }
 348     snew(atc->spline, atc->nthread);
 349     for (thread = 0; thread < atc->nthread; thread++)
 350     {
 351         if (atc->nthread > 1)
 352         {
 353             snew(atc->thread_plist[thread].n, atc->nthread+2*gmxCacheLineSize);
 354             atc->thread_plist[thread].n += gmxCacheLineSize;
 355         }
 356     }
 357 }
 358
 359 /*! \brief Destroy an atom communication data structure and its child structs */
 360 static void destroy_atomcomm(pme_atomcomm_t *atc)
 361 {
 362     sfree(atc->pd);
 363     if (atc->nslab > 1)
 364     {
 365         sfree(atc->node_dest);
 366         sfree(atc->node_src);
 367         for (int i = 0; i < atc->nthread; i++)
 368         {
 369             sfree(atc->count_thread[i]);
 370         }
 371         sfree(atc->count_thread);
 372         sfree(atc->rcount);
 373         sfree(atc->buf_index);
 374
 375         sfree(atc->x);
 376         sfree(atc->coefficient);
 377         sfree(atc->f);
 378     }
 379     sfree(atc->idx);
 380     sfree(atc->fractx);
 381
 382     sfree(atc->thread_idx);
 383     for (int i = 0; i < atc->nthread; i++)
 384     {
 385         if (atc->nthread > 1)
 386         {
 387             int *n_ptr = atc->thread_plist[i].n - gmxCacheLineSize;
 388             sfree(n_ptr);
 389             sfree(atc->thread_plist[i].i);
 390         }
 391         sfree(atc->spline[i].ind);
 392         for (int d = 0; d < ZZ; d++)
 393         {
 394             sfree(atc->spline[i].theta[d]);
 395             sfree(atc->spline[i].dtheta[d]);
 396         }
 397         sfree_aligned(atc->spline[i].ptr_dtheta_z);
 398         sfree_aligned(atc->spline[i].ptr_theta_z);
 399     }
 400     if (atc->nthread > 1)
 401     {
 402         sfree(atc->thread_plist);
 403     }
 404     sfree(atc->spline);
 405 }
 406
 407 /*! \brief Initialize data structure for communication */
 408 static void
 409 init_overlap_comm(pme_overlap_t *  ol,
 410                   int              norder,
 411 #if GMX_MPI
 412                   MPI_Comm         comm,
 413 #endif
 414                   int              nnodes,
 415                   int              nodeid,
 416                   int              ndata,
 417                   int              commplainsize)
 418 {
 419     gmx_bool         bCont;
 420 #if GMX_MPI
 421     MPI_Status       stat;
 422
 423     ol->mpi_comm = comm;
 424 #endif
 425
 426     ol->nnodes = nnodes;
 427     ol->nodeid = nodeid;
 428
 429     /* Linear translation of the PME grid won't affect reciprocal space
 430      * calculations, so to optimize we only interpolate "upwards",
 431      * which also means we only have to consider overlap in one direction.
 432      * I.e., particles on this node might also be spread to grid indices
 433      * that belong to higher nodes (modulo nnodes)
 434      */
 435
 436     ol->s2g0.resize(ol->nnodes + 1);
 437     ol->s2g1.resize(ol->nnodes);
 438     if (debug)
 439     {
 440         fprintf(debug, "PME slab boundaries:");
 441     }
 442     for (int i = 0; i < nnodes; i++)
 443     {
 444         /* s2g0 the local interpolation grid start.
 445          * s2g1 the local interpolation grid end.
 446          * Since in calc_pidx we divide particles, and not grid lines,
 447          * spatially uniform along dimension x or y, we need to round
 448          * s2g0 down and s2g1 up.
 449          */
 450         ol->s2g0[i] = (i * ndata + 0) / nnodes;
 451         ol->s2g1[i] = ((i + 1) * ndata + nnodes - 1) / nnodes + norder - 1;
 452
 453         if (debug)
 454         {
 455             fprintf(debug, "  %3d %3d", ol->s2g0[i], ol->s2g1[i]);
 456         }
 457     }
 458     ol->s2g0[nnodes] = ndata;
 459     if (debug)
 460     {
 461         fprintf(debug, "\n");
 462     }
 463
 464     /* Determine with how many nodes we need to communicate the grid overlap */
 465     int testRankCount = 0;
 466     do
 467     {
 468         testRankCount++;
 469         bCont = FALSE;
 470         for (int i = 0; i < nnodes; i++)
 471         {
 472             if ((i + testRankCount <  nnodes && ol->s2g1[i] > ol->s2g0[i + testRankCount]) ||
 473                 (i + testRankCount >= nnodes && ol->s2g1[i] > ol->s2g0[i + testRankCount - nnodes] + ndata))
 474             {
 475                 bCont = TRUE;
 476             }
 477         }
 478     }
 479     while (bCont && testRankCount < nnodes);
 480
 481     ol->comm_data.resize(testRankCount - 1);
 482     ol->send_size = 0;
 483
 484     for (size_t b = 0; b < ol->comm_data.size(); b++)
 485     {
 486         pme_grid_comm_t *pgc = &ol->comm_data[b];
 487
 488         /* Send */
 489         pgc->send_id = (ol->nodeid + (b + 1)) % ol->nnodes;
 490         int fft_start = ol->s2g0[pgc->send_id];
 491         int fft_end   = ol->s2g0[pgc->send_id + 1];
 492         if (pgc->send_id < nodeid)
 493         {
 494             fft_start += ndata;
 495             fft_end   += ndata;
 496         }
 497         int send_index1  = ol->s2g1[nodeid];
 498         send_index1      = std::min(send_index1, fft_end);
 499         pgc->send_index0 = fft_start;
 500         pgc->send_nindex = std::max(0, send_index1 - pgc->send_index0);
 501         ol->send_size   += pgc->send_nindex;
 502
 503         /* We always start receiving to the first index of our slab */
 504         pgc->recv_id     = (ol->nodeid - (b + 1) + ol->nnodes) % ol->nnodes;
 505         fft_start        = ol->s2g0[ol->nodeid];
 506         fft_end          = ol->s2g0[ol->nodeid + 1];
 507         int recv_index1  = ol->s2g1[pgc->recv_id];
 508         if (pgc->recv_id > nodeid)
 509         {
 510             recv_index1 -= ndata;
 511         }
 512         recv_index1      = std::min(recv_index1, fft_end);
 513         pgc->recv_index0 = fft_start;
 514         pgc->recv_nindex = std::max(0, recv_index1 - pgc->recv_index0);
 515     }
 516
 517 #if GMX_MPI
 518     /* Communicate the buffer sizes to receive */
 519     for (size_t b = 0; b < ol->comm_data.size(); b++)
 520     {
 521         MPI_Sendrecv(&ol->send_size, 1, MPI_INT, ol->comm_data[b].send_id, b,
 522                      &ol->comm_data[b].recv_size, 1, MPI_INT, ol->comm_data[b].recv_id, b,
 523                      ol->mpi_comm, &stat);
 524     }
 525 #endif
 526
 527     /* For non-divisible grid we need pme_order iso pme_order-1 */
 528     ol->sendbuf.resize(norder * commplainsize);
 529     ol->recvbuf.resize(norder * commplainsize);
 530 }
 531
 532 int minimalPmeGridSize(int pmeOrder)
 533 {
 534     /* The actual grid size limitations are:
 535      *   serial:        >= pme_order
 536      *   DD, no OpenMP: >= 2*(pme_order - 1)
 537      *   DD, OpenMP:    >= pme_order + 1
 538      * But we use the maximum for simplicity since in practice there is not
 539      * much performance difference between pme_order and 2*(pme_order -1).
 540      */
 541     int minimalSize = 2*(pmeOrder - 1);
 542
 543     GMX_RELEASE_ASSERT(pmeOrder >= 3, "pmeOrder has to be >= 3");
 544     GMX_RELEASE_ASSERT(minimalSize >= pmeOrder + 1, "The grid size should be >= pmeOrder + 1");
 545
 546     return minimalSize;
 547 }
 548
 549 bool gmx_pme_check_restrictions(int pme_order,
 550                                 int nkx, int nky, int nkz,
 551                                 int numPmeDomainsAlongX,
 552                                 bool useThreads,
 553                                 bool errorsAreFatal)
 554 {
 555     if (pme_order > PME_ORDER_MAX)
 556     {
 557         if (!errorsAreFatal)
 558         {
 559             return false;
 560         }
 561
 562         std::string message = gmx::formatString(
 563                     "pme_order (%d) is larger than the maximum allowed value (%d). Modify and recompile the code if you really need such a high order.",
 564                     pme_order, PME_ORDER_MAX);
 565         GMX_THROW(gmx::InconsistentInputError(message));
 566     }
 567
 568     const int minGridSize = minimalPmeGridSize(pme_order);
 569     if (nkx < minGridSize ||
 570         nky < minGridSize ||
 571         nkz < minGridSize)
 572     {
 573         if (!errorsAreFatal)
 574         {
 575             return false;
 576         }
 577         std::string message = gmx::formatString(
 578                     "The PME grid sizes need to be >= 2*(pme_order-1) (%d)",
 579                     minGridSize);
 580         GMX_THROW(gmx::InconsistentInputError(message));
 581     }
 582
 583     /* Check for a limitation of the (current) sum_fftgrid_dd code.
 584      * We only allow multiple communication pulses in dim 1, not in dim 0.
 585      */
 586     if (useThreads && (nkx < numPmeDomainsAlongX*pme_order &&
 587                        nkx != numPmeDomainsAlongX*(pme_order - 1)))
 588     {
 589         if (!errorsAreFatal)
 590         {
 591             return false;
 592         }
 593         gmx_fatal(FARGS, "The number of PME grid lines per rank along x is %g. But when using OpenMP threads, the number of grid lines per rank along x should be >= pme_order (%d) or = pmeorder-1. To resolve this issue, use fewer ranks along x (and possibly more along y and/or z) by specifying -dd manually.",
 594                   nkx/static_cast<double>(numPmeDomainsAlongX), pme_order);
 595     }
 596
 597     return true;
 598 }
 599
 600 /*! \brief Round \p enumerator */
 601 static int div_round_up(int enumerator, int denominator)
 602 {
 603     return (enumerator + denominator - 1)/denominator;
 604 }
 605
 606 gmx_pme_t *gmx_pme_init(const t_commrec         *cr,
 607                         const NumPmeDomains     &numPmeDomains,
 608                         const t_inputrec        *ir,
 609                         int                      homenr,
 610                         gmx_bool                 bFreeEnergy_q,
 611                         gmx_bool                 bFreeEnergy_lj,
 612                         gmx_bool                 bReproducible,
 613                         real                     ewaldcoeff_q,
 614                         real                     ewaldcoeff_lj,
 615                         int                      nthread,
 616                         PmeRunMode               runMode,
 617                         PmeGpu                  *pmeGpu,
 618                         const gmx_device_info_t *gpuInfo,
 619                         PmeGpuProgramHandle      pmeGpuProgram,
 620                         const gmx::MDLogger      & /*mdlog*/)
 621 {
 622     int               use_threads, sum_use_threads, i;
 623     ivec              ndata;
 624
 625     if (debug)
 626     {
 627         fprintf(debug, "Creating PME data structures.\n");
 628     }
 629
 630     gmx::unique_cptr<gmx_pme_t, gmx_pme_destroy> pme(new gmx_pme_t());
 631
 632     pme->sum_qgrid_tmp       = nullptr;
 633     pme->sum_qgrid_dd_tmp    = nullptr;
 634
 635     pme->buf_nalloc          = 0;
 636
 637     pme->nnodes              = 1;
 638     pme->bPPnode             = TRUE;
 639
 640     pme->nnodes_major        = numPmeDomains.x;
 641     pme->nnodes_minor        = numPmeDomains.y;
 642
 643 #if GMX_MPI
 644     if (numPmeDomains.x*numPmeDomains.y > 1)
 645     {
 646         pme->mpi_comm = cr->mpi_comm_mygroup;
 647
 648         MPI_Comm_rank(pme->mpi_comm, &pme->nodeid);
 649         MPI_Comm_size(pme->mpi_comm, &pme->nnodes);
 650         if (pme->nnodes != numPmeDomains.x*numPmeDomains.y)
 651         {
 652             gmx_incons("PME rank count mismatch");
 653         }
 654     }
 655     else
 656     {
 657         pme->mpi_comm = MPI_COMM_NULL;
 658     }
 659 #endif
 660
 661     if (pme->nnodes == 1)
 662     {
 663 #if GMX_MPI
 664         pme->mpi_comm_d[0] = MPI_COMM_NULL;
 665         pme->mpi_comm_d[1] = MPI_COMM_NULL;
 666 #endif
 667         pme->ndecompdim   = 0;
 668         pme->nodeid_major = 0;
 669         pme->nodeid_minor = 0;
 670 #if GMX_MPI
 671         pme->mpi_comm_d[0] = pme->mpi_comm_d[1] = MPI_COMM_NULL;
 672 #endif
 673     }
 674     else
 675     {
 676         if (numPmeDomains.y == 1)
 677         {
 678 #if GMX_MPI
 679             pme->mpi_comm_d[0] = pme->mpi_comm;
 680             pme->mpi_comm_d[1] = MPI_COMM_NULL;
 681 #endif
 682             pme->ndecompdim   = 1;
 683             pme->nodeid_major = pme->nodeid;
 684             pme->nodeid_minor = 0;
 685
 686         }
 687         else if (numPmeDomains.x == 1)
 688         {
 689 #if GMX_MPI
 690             pme->mpi_comm_d[0] = MPI_COMM_NULL;
 691             pme->mpi_comm_d[1] = pme->mpi_comm;
 692 #endif
 693             pme->ndecompdim   = 1;
 694             pme->nodeid_major = 0;
 695             pme->nodeid_minor = pme->nodeid;
 696         }
 697         else
 698         {
 699             if (pme->nnodes % numPmeDomains.x != 0)
 700             {
 701                 gmx_incons("For 2D PME decomposition, #PME ranks must be divisible by the number of domains along x");
 702             }
 703             pme->ndecompdim = 2;
 704
 705 #if GMX_MPI
 706             MPI_Comm_split(pme->mpi_comm, pme->nodeid % numPmeDomains.y,
 707                            pme->nodeid, &pme->mpi_comm_d[0]);  /* My communicator along major dimension */
 708             MPI_Comm_split(pme->mpi_comm, pme->nodeid/numPmeDomains.y,
 709                            pme->nodeid, &pme->mpi_comm_d[1]);  /* My communicator along minor dimension */
 710
 711             MPI_Comm_rank(pme->mpi_comm_d[0], &pme->nodeid_major);
 712             MPI_Comm_size(pme->mpi_comm_d[0], &pme->nnodes_major);
 713             MPI_Comm_rank(pme->mpi_comm_d[1], &pme->nodeid_minor);
 714             MPI_Comm_size(pme->mpi_comm_d[1], &pme->nnodes_minor);
 715 #endif
 716         }
 717         pme->bPPnode = thisRankHasDuty(cr, DUTY_PP);
 718     }
 719
 720     pme->nthread = nthread;
 721
 722     /* Check if any of the PME MPI ranks uses threads */
 723     use_threads = (pme->nthread > 1 ? 1 : 0);
 724 #if GMX_MPI
 725     if (pme->nnodes > 1)
 726     {
 727         MPI_Allreduce(&use_threads, &sum_use_threads, 1, MPI_INT,
 728                       MPI_SUM, pme->mpi_comm);
 729     }
 730     else
 731 #endif
 732     {
 733         sum_use_threads = use_threads;
 734     }
 735     pme->bUseThreads = (sum_use_threads > 0);
 736
 737     if (ir->ePBC == epbcSCREW)
 738     {
 739         gmx_fatal(FARGS, "pme does not (yet) work with pbc = screw");
 740     }
 741
 742     /* NOTE:
 743      * It is likely that the current gmx_pme_do() routine supports calculating
 744      * only Coulomb or LJ while gmx_pme_init() configures for both,
 745      * but that has never been tested.
 746      * It is likely that the current gmx_pme_do() routine supports calculating,
 747      * not calculating free-energy for Coulomb and/or LJ while gmx_pme_init()
 748      * configures with free-energy, but that has never been tested.
 749      */
 750     pme->doCoulomb     = EEL_PME(ir->coulombtype);
 751     pme->doLJ          = EVDW_PME(ir->vdwtype);
 752     pme->bFEP_q        = ((ir->efep != efepNO) && bFreeEnergy_q);
 753     pme->bFEP_lj       = ((ir->efep != efepNO) && bFreeEnergy_lj);
 754     pme->bFEP          = (pme->bFEP_q || pme->bFEP_lj);
 755     pme->nkx           = ir->nkx;
 756     pme->nky           = ir->nky;
 757     pme->nkz           = ir->nkz;
 758     pme->bP3M          = (ir->coulombtype == eelP3M_AD || getenv("GMX_PME_P3M") != nullptr);
 759     pme->pme_order     = ir->pme_order;
 760     pme->ewaldcoeff_q  = ewaldcoeff_q;
 761     pme->ewaldcoeff_lj = ewaldcoeff_lj;
 762
 763     /* Always constant electrostatics coefficients */
 764     pme->epsilon_r     = ir->epsilon_r;
 765
 766     /* Always constant LJ coefficients */
 767     pme->ljpme_combination_rule = ir->ljpme_combination_rule;
 768
 769     // The box requires scaling with nwalls = 2, we store that condition as well
 770     // as the scaling factor
 771     delete pme->boxScaler;
 772     pme->boxScaler = new EwaldBoxZScaler(*ir);
 773
 774     /* If we violate restrictions, generate a fatal error here */
 775     gmx_pme_check_restrictions(pme->pme_order,
 776                                pme->nkx, pme->nky, pme->nkz,
 777                                pme->nnodes_major,
 778                                pme->bUseThreads,
 779                                true);
 780
 781     if (pme->nnodes > 1)
 782     {
 783         double imbal;
 784
 785 #if GMX_MPI
 786         MPI_Type_contiguous(DIM, GMX_MPI_REAL, &(pme->rvec_mpi));
 787         MPI_Type_commit(&(pme->rvec_mpi));
 788 #endif
 789
 790         /* Note that the coefficient spreading and force gathering, which usually
 791          * takes about the same amount of time as FFT+solve_pme,
 792          * is always fully load balanced
 793          * (unless the coefficient distribution is inhomogeneous).
 794          */
 795
 796         imbal = estimate_pme_load_imbalance(pme.get());
 797         if (imbal >= 1.2 && pme->nodeid_major == 0 && pme->nodeid_minor == 0)
 798         {
 799             fprintf(stderr,
 800                     "\n"
 801                     "NOTE: The load imbalance in PME FFT and solve is %d%%.\n"
 802                     "      For optimal PME load balancing\n"
 803                     "      PME grid_x (%d) and grid_y (%d) should be divisible by #PME_ranks_x (%d)\n"
 804                     "      and PME grid_y (%d) and grid_z (%d) should be divisible by #PME_ranks_y (%d)\n"
 805                     "\n",
 806                     gmx::roundToInt((imbal-1)*100),
 807                     pme->nkx, pme->nky, pme->nnodes_major,
 808                     pme->nky, pme->nkz, pme->nnodes_minor);
 809         }
 810     }
 811
 812     /* For non-divisible grid we need pme_order iso pme_order-1 */
 813     /* In sum_qgrid_dd x overlap is copied in place: take padding into account.
 814      * y is always copied through a buffer: we don't need padding in z,
 815      * but we do need the overlap in x because of the communication order.
 816      */
 817     init_overlap_comm(&pme->overlap[0], pme->pme_order,
 818 #if GMX_MPI
 819                       pme->mpi_comm_d[0],
 820 #endif
 821                       pme->nnodes_major, pme->nodeid_major,
 822                       pme->nkx,
 823                       (div_round_up(pme->nky, pme->nnodes_minor)+pme->pme_order)*(pme->nkz+pme->pme_order-1));
 824
 825     /* Along overlap dim 1 we can send in multiple pulses in sum_fftgrid_dd.
 826      * We do this with an offset buffer of equal size, so we need to allocate
 827      * extra for the offset. That's what the (+1)*pme->nkz is for.
 828      */
 829     init_overlap_comm(&pme->overlap[1], pme->pme_order,
 830 #if GMX_MPI
 831                       pme->mpi_comm_d[1],
 832 #endif
 833                       pme->nnodes_minor, pme->nodeid_minor,
 834                       pme->nky,
 835                       (div_round_up(pme->nkx, pme->nnodes_major)+pme->pme_order+1)*pme->nkz);
 836
 837     /* Double-check for a limitation of the (current) sum_fftgrid_dd code.
 838      * Note that gmx_pme_check_restrictions checked for this already.
 839      */
 840     if (pme->bUseThreads && (pme->overlap[0].comm_data.size() > 1))
 841     {
 842         gmx_incons("More than one communication pulse required for grid overlap communication along the major dimension while using threads");
 843     }
 844
 845     snew(pme->bsp_mod[XX], pme->nkx);
 846     snew(pme->bsp_mod[YY], pme->nky);
 847     snew(pme->bsp_mod[ZZ], pme->nkz);
 848
 849     pme->gpu     = pmeGpu; /* Carrying over the single GPU structure */
 850     pme->runMode = runMode;
 851
 852     /* The required size of the interpolation grid, including overlap.
 853      * The allocated size (pmegrid_n?) might be slightly larger.
 854      */
 855     pme->pmegrid_nx = pme->overlap[0].s2g1[pme->nodeid_major] -
 856         pme->overlap[0].s2g0[pme->nodeid_major];
 857     pme->pmegrid_ny = pme->overlap[1].s2g1[pme->nodeid_minor] -
 858         pme->overlap[1].s2g0[pme->nodeid_minor];
 859     pme->pmegrid_nz_base = pme->nkz;
 860     pme->pmegrid_nz      = pme->pmegrid_nz_base + pme->pme_order - 1;
 861     set_grid_alignment(&pme->pmegrid_nz, pme->pme_order);
 862     pme->pmegrid_start_ix = pme->overlap[0].s2g0[pme->nodeid_major];
 863     pme->pmegrid_start_iy = pme->overlap[1].s2g0[pme->nodeid_minor];
 864     pme->pmegrid_start_iz = 0;
 865
 866     make_gridindex_to_localindex(pme->nkx,
 867                                  pme->pmegrid_start_ix,
 868                                  pme->pmegrid_nx - (pme->pme_order-1),
 869                                  &pme->nnx, &pme->fshx);
 870     make_gridindex_to_localindex(pme->nky,
 871                                  pme->pmegrid_start_iy,
 872                                  pme->pmegrid_ny - (pme->pme_order-1),
 873                                  &pme->nny, &pme->fshy);
 874     make_gridindex_to_localindex(pme->nkz,
 875                                  pme->pmegrid_start_iz,
 876                                  pme->pmegrid_nz_base,
 877                                  &pme->nnz, &pme->fshz);
 878
 879     pme->spline_work = make_pme_spline_work(pme->pme_order);
 880
 881     ndata[0]    = pme->nkx;
 882     ndata[1]    = pme->nky;
 883     ndata[2]    = pme->nkz;
 884     /* It doesn't matter if we allocate too many grids here,
 885      * we only allocate and use the ones we need.
 886      */
 887     if (pme->doLJ)
 888     {
 889         pme->ngrids = ((ir->ljpme_combination_rule == eljpmeLB) ? DO_Q_AND_LJ_LB : DO_Q_AND_LJ);
 890     }
 891     else
 892     {
 893         pme->ngrids = DO_Q;
 894     }
 895     snew(pme->fftgrid, pme->ngrids);
 896     snew(pme->cfftgrid, pme->ngrids);
 897     snew(pme->pfft_setup, pme->ngrids);
 898
 899     for (i = 0; i < pme->ngrids; ++i)
 900     {
 901         if ((i <  DO_Q && pme->doCoulomb && (i == 0 ||
 902                                              bFreeEnergy_q)) ||
 903             (i >= DO_Q && pme->doLJ && (i == 2 ||
 904                                         bFreeEnergy_lj ||
 905                                         ir->ljpme_combination_rule == eljpmeLB)))
 906         {
 907             pmegrids_init(&pme->pmegrid[i],
 908                           pme->pmegrid_nx, pme->pmegrid_ny, pme->pmegrid_nz,
 909                           pme->pmegrid_nz_base,
 910                           pme->pme_order,
 911                           pme->bUseThreads,
 912                           pme->nthread,
 913                           pme->overlap[0].s2g1[pme->nodeid_major]-pme->overlap[0].s2g0[pme->nodeid_major+1],
 914                           pme->overlap[1].s2g1[pme->nodeid_minor]-pme->overlap[1].s2g0[pme->nodeid_minor+1]);
 915             /* This routine will allocate the grid data to fit the FFTs */
 916             const auto allocateRealGridForGpu = (pme->runMode == PmeRunMode::Mixed) ? gmx::PinningPolicy::PinnedIfSupported : gmx::PinningPolicy::CannotBePinned;
 917             gmx_parallel_3dfft_init(&pme->pfft_setup[i], ndata,
 918                                     &pme->fftgrid[i], &pme->cfftgrid[i],
 919                                     pme->mpi_comm_d,
 920                                     bReproducible, pme->nthread, allocateRealGridForGpu);
 921
 922         }
 923     }
 924
 925     if (!pme->bP3M)
 926     {
 927         /* Use plain SPME B-spline interpolation */
 928         make_bspline_moduli(pme->bsp_mod, pme->nkx, pme->nky, pme->nkz, pme->pme_order);
 929     }
 930     else
 931     {
 932         /* Use the P3M grid-optimized influence function */
 933         make_p3m_bspline_moduli(pme->bsp_mod, pme->nkx, pme->nky, pme->nkz, pme->pme_order);
 934     }
 935
 936     /* Use atc[0] for spreading */
 937     init_atomcomm(pme.get(), &pme->atc[0], numPmeDomains.x > 1 ? 0 : 1, TRUE);
 938     if (pme->ndecompdim >= 2)
 939     {
 940         init_atomcomm(pme.get(), &pme->atc[1], 1, FALSE);
 941     }
 942
 943     if (pme->nnodes == 1)
 944     {
 945         pme->atc[0].n = homenr;
 946         pme_realloc_atomcomm_things(&pme->atc[0]);
 947     }
 948
 949     pme->lb_buf1       = nullptr;
 950     pme->lb_buf2       = nullptr;
 951     pme->lb_buf_nalloc = 0;
 952
 953     if (pme_gpu_active(pme.get()))
 954     {
 955         if (!pme->gpu)
 956         {
 957             // Initial check of validity of the data
 958             std::string errorString;
 959             bool        canRunOnGpu = pme_gpu_check_restrictions(pme.get(), &errorString);
 960             if (!canRunOnGpu)
 961             {
 962                 GMX_THROW(gmx::NotImplementedError(errorString));
 963             }
 964         }
 965
 966         pme_gpu_reinit(pme.get(), gpuInfo, pmeGpuProgram);
 967     }
 968
 969     pme_init_all_work(&pme->solve_work, pme->nthread, pme->nkx);
 970
 971     // no exception was thrown during the init, so we hand over the PME structure handle
 972     return pme.release();
 973 }
 974
 975 void gmx_pme_reinit(struct gmx_pme_t **pmedata,
 976                     const t_commrec   *cr,
 977                     struct gmx_pme_t * pme_src,
 978                     const t_inputrec * ir,
 979                     const ivec         grid_size,
 980                     real               ewaldcoeff_q,
 981                     real               ewaldcoeff_lj)
 982 {
 983     int        homenr;
 984
 985     // Create a copy of t_inputrec fields that are used in gmx_pme_init().
 986     // TODO: This would be better as just copying a sub-structure that contains
 987     // all the PME parameters and nothing else.
 988     t_inputrec irc;
 989     irc.ePBC                   = ir->ePBC;
 990     irc.coulombtype            = ir->coulombtype;
 991     irc.vdwtype                = ir->vdwtype;
 992     irc.efep                   = ir->efep;
 993     irc.pme_order              = ir->pme_order;
 994     irc.epsilon_r              = ir->epsilon_r;
 995     irc.ljpme_combination_rule = ir->ljpme_combination_rule;
 996     irc.nkx                    = grid_size[XX];
 997     irc.nky                    = grid_size[YY];
 998     irc.nkz                    = grid_size[ZZ];
 999
1000     if (pme_src->nnodes == 1)
1001     {
1002         homenr = pme_src->atc[0].n;
1003     }
1004     else
1005     {
1006         homenr = -1;
1007     }
1008
1009     try
1010     {
1011         const gmx::MDLogger dummyLogger;
1012         // This is reinit which is currently only changing grid size/coefficients,
1013         // so we don't expect the actual logging.
1014         // TODO: when PME is an object, it should take reference to mdlog on construction and save it.
1015         GMX_ASSERT(pmedata, "Invalid PME pointer");
1016         NumPmeDomains numPmeDomains = { pme_src->nnodes_major, pme_src->nnodes_minor };
1017         *pmedata = gmx_pme_init(cr, numPmeDomains,
1018                                 &irc, homenr, pme_src->bFEP_q, pme_src->bFEP_lj, FALSE, ewaldcoeff_q, ewaldcoeff_lj,
1019                                 pme_src->nthread, pme_src->runMode, pme_src->gpu, nullptr, nullptr, dummyLogger);
1020         //TODO this is mostly passing around current values
1021     }
1022     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1023
1024     /* We can easily reuse the allocated pme grids in pme_src */
1025     reuse_pmegrids(&pme_src->pmegrid[PME_GRID_QA], &(*pmedata)->pmegrid[PME_GRID_QA]);
1026     /* We would like to reuse the fft grids, but that's harder */
1027 }
1028
1029 void gmx_pme_calc_energy(struct gmx_pme_t *pme, int n, rvec *x, real *q, real *V)
1030 {
1031     pme_atomcomm_t *atc;
1032     pmegrids_t     *grid;
1033
1034     if (pme->nnodes > 1)
1035     {
1036         gmx_incons("gmx_pme_calc_energy called in parallel");
1037     }
1038     if (pme->bFEP_q)
1039     {
1040         gmx_incons("gmx_pme_calc_energy with free energy");
1041     }
1042
1043     atc            = &pme->atc_energy;
1044     atc->nthread   = 1;
1045     if (atc->spline == nullptr)
1046     {
1047         snew(atc->spline, atc->nthread);
1048     }
1049     atc->nslab     = 1;
1050     atc->bSpread   = TRUE;
1051     atc->pme_order = pme->pme_order;
1052     atc->n         = n;
1053     pme_realloc_atomcomm_things(atc);
1054     atc->x           = x;
1055     atc->coefficient = q;
1056
1057     /* We only use the A-charges grid */
1058     grid = &pme->pmegrid[PME_GRID_QA];
1059
1060     /* Only calculate the spline coefficients, don't actually spread */
1061     spread_on_grid(pme, atc, nullptr, TRUE, FALSE, pme->fftgrid[PME_GRID_QA], FALSE, PME_GRID_QA);
1062
1063     *V = gather_energy_bsplines(pme, grid->grid.grid, atc);
1064 }
1065
1066 /*! \brief Calculate initial Lorentz-Berthelot coefficients for LJ-PME */
1067 static void
1068 calc_initial_lb_coeffs(struct gmx_pme_t *pme, const real *local_c6, const real *local_sigma)
1069 {
1070     int  i;
1071     for (i = 0; i < pme->atc[0].n; ++i)
1072     {
1073         real sigma4;
1074         sigma4                     = local_sigma[i];
1075         sigma4                     = sigma4*sigma4;
1076         sigma4                     = sigma4*sigma4;
1077         pme->atc[0].coefficient[i] = local_c6[i] / sigma4;
1078     }
1079 }
1080
1081 /*! \brief Calculate next Lorentz-Berthelot coefficients for LJ-PME */
1082 static void
1083 calc_next_lb_coeffs(struct gmx_pme_t *pme, const real *local_sigma)
1084 {
1085     int  i;
1086
1087     for (i = 0; i < pme->atc[0].n; ++i)
1088     {
1089         pme->atc[0].coefficient[i] *= local_sigma[i];
1090     }
1091 }
1092
1093 int gmx_pme_do(struct gmx_pme_t *pme,
1094                int start,       int homenr,
1095                rvec x[],        rvec f[],
1096                real chargeA[],  real chargeB[],
1097                real c6A[],      real c6B[],
1098                real sigmaA[],   real sigmaB[],
1099                matrix box,      const t_commrec *cr,
1100                int  maxshift_x, int maxshift_y,
1101                t_nrnb *nrnb,    gmx_wallcycle *wcycle,
1102                matrix vir_q,    matrix vir_lj,
1103                real *energy_q,  real *energy_lj,
1104                real lambda_q,   real lambda_lj,
1105                real *dvdlambda_q, real *dvdlambda_lj,
1106                int flags)
1107 {
1108     GMX_ASSERT(pme->runMode == PmeRunMode::CPU, "gmx_pme_do should not be called on the GPU PME run.");
1109
1110     int                  d, i, j, npme, grid_index, max_grid_index;
1111     int                  n_d;
1112     pme_atomcomm_t      *atc        = nullptr;
1113     pmegrids_t          *pmegrid    = nullptr;
1114     real                *grid       = nullptr;
1115     rvec                *f_d;
1116     real                *coefficient = nullptr;
1117     real                 energy_AB[4];
1118     matrix               vir_AB[4];
1119     real                 scale, lambda;
1120     gmx_bool             bClearF;
1121     gmx_parallel_3dfft_t pfft_setup;
1122     real              *  fftgrid;
1123     t_complex          * cfftgrid;
1124     int                  thread;
1125     gmx_bool             bFirst, bDoSplines;
1126     int                  fep_state;
1127     int                  fep_states_lj           = pme->bFEP_lj ? 2 : 1;
1128     const gmx_bool       bCalcEnerVir            = (flags & GMX_PME_CALC_ENER_VIR) != 0;
1129     const gmx_bool       bBackFFT                = (flags & (GMX_PME_CALC_F | GMX_PME_CALC_POT)) != 0;
1130     const gmx_bool       bCalcF                  = (flags & GMX_PME_CALC_F) != 0;
1131
1132     /* We could be passing lambda!=1 while no q or LJ is actually perturbed */
1133     if (!pme->bFEP_q)
1134     {
1135         lambda_q  = 1;
1136     }
1137     if (!pme->bFEP_lj)
1138     {
1139         lambda_lj = 1;
1140     }
1141
1142     assert(pme->nnodes > 0);
1143     assert(pme->nnodes == 1 || pme->ndecompdim > 0);
1144
1145     if (pme->nnodes > 1)
1146     {
1147         atc      = &pme->atc[0];
1148         atc->npd = homenr;
1149         if (atc->npd > atc->pd_nalloc)
1150         {
1151             atc->pd_nalloc = over_alloc_dd(atc->npd);
1152             srenew(atc->pd, atc->pd_nalloc);
1153         }
1154         for (d = pme->ndecompdim-1; d >= 0; d--)
1155         {
1156             atc           = &pme->atc[d];
1157             atc->maxshift = (atc->dimind == 0 ? maxshift_x : maxshift_y);
1158         }
1159     }
1160     else
1161     {
1162         atc = &pme->atc[0];
1163         /* This could be necessary for TPI */
1164         pme->atc[0].n = homenr;
1165         if (DOMAINDECOMP(cr))
1166         {
1167             pme_realloc_atomcomm_things(atc);
1168         }
1169         atc->x = x;
1170         atc->f = f;
1171     }
1172
1173     matrix scaledBox;
1174     pme->boxScaler->scaleBox(box, scaledBox);
1175
1176     gmx::invertBoxMatrix(scaledBox, pme->recipbox);
1177     bFirst = TRUE;
1178
1179     /* For simplicity, we construct the splines for all particles if
1180      * more than one PME calculations is needed. Some optimization
1181      * could be done by keeping track of which atoms have splines
1182      * constructed, and construct new splines on each pass for atoms
1183      * that don't yet have them.
1184      */
1185
1186     bDoSplines = pme->bFEP || (pme->doCoulomb && pme->doLJ);
1187
1188     /* We need a maximum of four separate PME calculations:
1189      * grid_index=0: Coulomb PME with charges from state A
1190      * grid_index=1: Coulomb PME with charges from state B
1191      * grid_index=2: LJ PME with C6 from state A
1192      * grid_index=3: LJ PME with C6 from state B
1193      * For Lorentz-Berthelot combination rules, a separate loop is used to
1194      * calculate all the terms
1195      */
1196
1197     /* If we are doing LJ-PME with LB, we only do Q here */
1198     max_grid_index = (pme->ljpme_combination_rule == eljpmeLB) ? DO_Q : DO_Q_AND_LJ;
1199
1200     for (grid_index = 0; grid_index < max_grid_index; ++grid_index)
1201     {
1202         /* Check if we should do calculations at this grid_index
1203          * If grid_index is odd we should be doing FEP
1204          * If grid_index < 2 we should be doing electrostatic PME
1205          * If grid_index >= 2 we should be doing LJ-PME
1206          */
1207         if ((grid_index <  DO_Q && (!pme->doCoulomb ||
1208                                     (grid_index == 1 && !pme->bFEP_q))) ||
1209             (grid_index >= DO_Q && (!pme->doLJ ||
1210                                     (grid_index == 3 && !pme->bFEP_lj))))
1211         {
1212             continue;
1213         }
1214         /* Unpack structure */
1215         pmegrid    = &pme->pmegrid[grid_index];
1216         fftgrid    = pme->fftgrid[grid_index];
1217         cfftgrid   = pme->cfftgrid[grid_index];
1218         pfft_setup = pme->pfft_setup[grid_index];
1219         switch (grid_index)
1220         {
1221             case 0: coefficient = chargeA + start; break;
1222             case 1: coefficient = chargeB + start; break;
1223             case 2: coefficient = c6A + start; break;
1224             case 3: coefficient = c6B + start; break;
1225         }
1226
1227         grid = pmegrid->grid.grid;
1228
1229         if (debug)
1230         {
1231             fprintf(debug, "PME: number of ranks = %d, rank = %d\n",
1232                     cr->nnodes, cr->nodeid);
1233             fprintf(debug, "Grid = %p\n", static_cast<void*>(grid));
1234             if (grid == nullptr)
1235             {
1236                 gmx_fatal(FARGS, "No grid!");
1237             }
1238         }
1239
1240         if (pme->nnodes == 1)
1241         {
1242             atc->coefficient = coefficient;
1243         }
1244         else
1245         {
1246             wallcycle_start(wcycle, ewcPME_REDISTXF);
1247             do_redist_pos_coeffs(pme, cr, start, homenr, bFirst, x, coefficient);
1248
1249             wallcycle_stop(wcycle, ewcPME_REDISTXF);
1250         }
1251
1252         if (debug)
1253         {
1254             fprintf(debug, "Rank= %6d, pme local particles=%6d\n",
1255                     cr->nodeid, atc->n);
1256         }
1257
1258         if (flags & GMX_PME_SPREAD)
1259         {
1260             wallcycle_start(wcycle, ewcPME_SPREAD);
1261
1262             /* Spread the coefficients on a grid */
1263             spread_on_grid(pme, &pme->atc[0], pmegrid, bFirst, TRUE, fftgrid, bDoSplines, grid_index);
1264
1265             if (bFirst)
1266             {
1267                 inc_nrnb(nrnb, eNR_WEIGHTS, DIM*atc->n);
1268             }
1269             inc_nrnb(nrnb, eNR_SPREADBSP,
1270                      pme->pme_order*pme->pme_order*pme->pme_order*atc->n);
1271
1272             if (!pme->bUseThreads)
1273             {
1274                 wrap_periodic_pmegrid(pme, grid);
1275
1276                 /* sum contributions to local grid from other nodes */
1277 #if GMX_MPI
1278                 if (pme->nnodes > 1)
1279                 {
1280                     gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_FORWARD);
1281                 }
1282 #endif
1283
1284                 copy_pmegrid_to_fftgrid(pme, grid, fftgrid, grid_index);
1285             }
1286
1287             wallcycle_stop(wcycle, ewcPME_SPREAD);
1288
1289             /* TODO If the OpenMP and single-threaded implementations
1290                converge, then spread_on_grid() and
1291                copy_pmegrid_to_fftgrid() will perhaps live in the same
1292                source file.
1293              */
1294         }
1295
1296         /* Here we start a large thread parallel region */
1297 #pragma omp parallel num_threads(pme->nthread) private(thread)
1298         {
1299             try
1300             {
1301                 thread = gmx_omp_get_thread_num();
1302                 if (flags & GMX_PME_SOLVE)
1303                 {
1304                     int loop_count;
1305
1306                     /* do 3d-fft */
1307                     if (thread == 0)
1308                     {
1309                         wallcycle_start(wcycle, ewcPME_FFT);
1310                     }
1311                     gmx_parallel_3dfft_execute(pfft_setup, GMX_FFT_REAL_TO_COMPLEX,
1312                                                thread, wcycle);
1313                     if (thread == 0)
1314                     {
1315                         wallcycle_stop(wcycle, ewcPME_FFT);
1316                     }
1317
1318                     /* solve in k-space for our local cells */
1319                     if (thread == 0)
1320                     {
1321                         wallcycle_start(wcycle, (grid_index < DO_Q ? ewcPME_SOLVE : ewcLJPME));
1322                     }
1323                     if (grid_index < DO_Q)
1324                     {
1325                         loop_count =
1326                             solve_pme_yzx(pme, cfftgrid,
1327                                           scaledBox[XX][XX]*scaledBox[YY][YY]*scaledBox[ZZ][ZZ],
1328                                           bCalcEnerVir,
1329                                           pme->nthread, thread);
1330                     }
1331                     else
1332                     {
1333                         loop_count =
1334                             solve_pme_lj_yzx(pme, &cfftgrid, FALSE,
1335                                              scaledBox[XX][XX]*scaledBox[YY][YY]*scaledBox[ZZ][ZZ],
1336                                              bCalcEnerVir,
1337                                              pme->nthread, thread);
1338                     }
1339
1340                     if (thread == 0)
1341                     {
1342                         wallcycle_stop(wcycle, (grid_index < DO_Q ? ewcPME_SOLVE : ewcLJPME));
1343                         inc_nrnb(nrnb, eNR_SOLVEPME, loop_count);
1344                     }
1345                 }
1346
1347                 if (bBackFFT)
1348                 {
1349                     /* do 3d-invfft */
1350                     if (thread == 0)
1351                     {
1352                         wallcycle_start(wcycle, ewcPME_FFT);
1353                     }
1354                     gmx_parallel_3dfft_execute(pfft_setup, GMX_FFT_COMPLEX_TO_REAL,
1355                                                thread, wcycle);
1356                     if (thread == 0)
1357                     {
1358                         wallcycle_stop(wcycle, ewcPME_FFT);
1359
1360
1361                         if (pme->nodeid == 0)
1362                         {
1363                             real ntot = pme->nkx*pme->nky*pme->nkz;
1364                             npme  = static_cast<int>(ntot*std::log(ntot)/std::log(2.0));
1365                             inc_nrnb(nrnb, eNR_FFT, 2*npme);
1366                         }
1367
1368                         /* Note: this wallcycle region is closed below
1369                            outside an OpenMP region, so take care if
1370                            refactoring code here. */
1371                         wallcycle_start(wcycle, ewcPME_GATHER);
1372                     }
1373
1374                     copy_fftgrid_to_pmegrid(pme, fftgrid, grid, grid_index, pme->nthread, thread);
1375                 }
1376             } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1377         }
1378         /* End of thread parallel section.
1379          * With MPI we have to synchronize here before gmx_sum_qgrid_dd.
1380          */
1381
1382         if (bBackFFT)
1383         {
1384             /* distribute local grid to all nodes */
1385 #if GMX_MPI
1386             if (pme->nnodes > 1)
1387             {
1388                 gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_BACKWARD);
1389             }
1390 #endif
1391
1392             unwrap_periodic_pmegrid(pme, grid);
1393         }
1394
1395         if (bCalcF)
1396         {
1397             /* interpolate forces for our local atoms */
1398
1399
1400             /* If we are running without parallelization,
1401              * atc->f is the actual force array, not a buffer,
1402              * therefore we should not clear it.
1403              */
1404             lambda  = grid_index < DO_Q ? lambda_q : lambda_lj;
1405             bClearF = (bFirst && PAR(cr));
1406 #pragma omp parallel for num_threads(pme->nthread) schedule(static)
1407             for (thread = 0; thread < pme->nthread; thread++)
1408             {
1409                 try
1410                 {
1411                     gather_f_bsplines(pme, grid, bClearF, atc,
1412                                       &atc->spline[thread],
1413                                       pme->bFEP ? (grid_index % 2 == 0 ? 1.0-lambda : lambda) : 1.0);
1414                 }
1415                 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1416             }
1417
1418
1419             inc_nrnb(nrnb, eNR_GATHERFBSP,
1420                      pme->pme_order*pme->pme_order*pme->pme_order*pme->atc[0].n);
1421             /* Note: this wallcycle region is opened above inside an OpenMP
1422                region, so take care if refactoring code here. */
1423             wallcycle_stop(wcycle, ewcPME_GATHER);
1424         }
1425
1426         if (bCalcEnerVir)
1427         {
1428             /* This should only be called on the master thread
1429              * and after the threads have synchronized.
1430              */
1431             if (grid_index < 2)
1432             {
1433                 get_pme_ener_vir_q(pme->solve_work, pme->nthread, &energy_AB[grid_index], vir_AB[grid_index]);
1434             }
1435             else
1436             {
1437                 get_pme_ener_vir_lj(pme->solve_work, pme->nthread, &energy_AB[grid_index], vir_AB[grid_index]);
1438             }
1439         }
1440         bFirst = FALSE;
1441     } /* of grid_index-loop */
1442
1443     /* For Lorentz-Berthelot combination rules in LJ-PME, we need to calculate
1444      * seven terms. */
1445
1446     if (pme->doLJ && pme->ljpme_combination_rule == eljpmeLB)
1447     {
1448         /* Loop over A- and B-state if we are doing FEP */
1449         for (fep_state = 0; fep_state < fep_states_lj; ++fep_state)
1450         {
1451             real *local_c6 = nullptr, *local_sigma = nullptr, *RedistC6 = nullptr, *RedistSigma = nullptr;
1452             if (pme->nnodes == 1)
1453             {
1454                 if (pme->lb_buf1 == nullptr)
1455                 {
1456                     pme->lb_buf_nalloc = pme->atc[0].n;
1457                     snew(pme->lb_buf1, pme->lb_buf_nalloc);
1458                 }
1459                 pme->atc[0].coefficient = pme->lb_buf1;
1460                 switch (fep_state)
1461                 {
1462                     case 0:
1463                         local_c6      = c6A;
1464                         local_sigma   = sigmaA;
1465                         break;
1466                     case 1:
1467                         local_c6      = c6B;
1468                         local_sigma   = sigmaB;
1469                         break;
1470                     default:
1471                         gmx_incons("Trying to access wrong FEP-state in LJ-PME routine");
1472                 }
1473             }
1474             else
1475             {
1476                 atc = &pme->atc[0];
1477                 switch (fep_state)
1478                 {
1479                     case 0:
1480                         RedistC6      = c6A;
1481                         RedistSigma   = sigmaA;
1482                         break;
1483                     case 1:
1484                         RedistC6      = c6B;
1485                         RedistSigma   = sigmaB;
1486                         break;
1487                     default:
1488                         gmx_incons("Trying to access wrong FEP-state in LJ-PME routine");
1489                 }
1490                 wallcycle_start(wcycle, ewcPME_REDISTXF);
1491
1492                 do_redist_pos_coeffs(pme, cr, start, homenr, bFirst, x, RedistC6);
1493                 if (pme->lb_buf_nalloc < atc->n)
1494                 {
1495                     pme->lb_buf_nalloc = atc->nalloc;
1496                     srenew(pme->lb_buf1, pme->lb_buf_nalloc);
1497                     srenew(pme->lb_buf2, pme->lb_buf_nalloc);
1498                 }
1499                 local_c6 = pme->lb_buf1;
1500                 for (i = 0; i < atc->n; ++i)
1501                 {
1502                     local_c6[i] = atc->coefficient[i];
1503                 }
1504
1505                 do_redist_pos_coeffs(pme, cr, start, homenr, FALSE, x, RedistSigma);
1506                 local_sigma = pme->lb_buf2;
1507                 for (i = 0; i < atc->n; ++i)
1508                 {
1509                     local_sigma[i] = atc->coefficient[i];
1510                 }
1511
1512                 wallcycle_stop(wcycle, ewcPME_REDISTXF);
1513             }
1514             calc_initial_lb_coeffs(pme, local_c6, local_sigma);
1515
1516             /*Seven terms in LJ-PME with LB, grid_index < 2 reserved for electrostatics*/
1517             for (grid_index = 2; grid_index < 9; ++grid_index)
1518             {
1519                 /* Unpack structure */
1520                 pmegrid    = &pme->pmegrid[grid_index];
1521                 fftgrid    = pme->fftgrid[grid_index];
1522                 pfft_setup = pme->pfft_setup[grid_index];
1523                 calc_next_lb_coeffs(pme, local_sigma);
1524                 grid = pmegrid->grid.grid;
1525
1526                 if (flags & GMX_PME_SPREAD)
1527                 {
1528                     wallcycle_start(wcycle, ewcPME_SPREAD);
1529                     /* Spread the c6 on a grid */
1530                     spread_on_grid(pme, &pme->atc[0], pmegrid, bFirst, TRUE, fftgrid, bDoSplines, grid_index);
1531
1532                     if (bFirst)
1533                     {
1534                         inc_nrnb(nrnb, eNR_WEIGHTS, DIM*atc->n);
1535                     }
1536
1537                     inc_nrnb(nrnb, eNR_SPREADBSP,
1538                              pme->pme_order*pme->pme_order*pme->pme_order*atc->n);
1539                     if (pme->nthread == 1)
1540                     {
1541                         wrap_periodic_pmegrid(pme, grid);
1542                         /* sum contributions to local grid from other nodes */
1543 #if GMX_MPI
1544                         if (pme->nnodes > 1)
1545                         {
1546                             gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_FORWARD);
1547                         }
1548 #endif
1549                         copy_pmegrid_to_fftgrid(pme, grid, fftgrid, grid_index);
1550                     }
1551                     wallcycle_stop(wcycle, ewcPME_SPREAD);
1552                 }
1553                 /*Here we start a large thread parallel region*/
1554 #pragma omp parallel num_threads(pme->nthread) private(thread)
1555                 {
1556                     try
1557                     {
1558                         thread = gmx_omp_get_thread_num();
1559                         if (flags & GMX_PME_SOLVE)
1560                         {
1561                             /* do 3d-fft */
1562                             if (thread == 0)
1563                             {
1564                                 wallcycle_start(wcycle, ewcPME_FFT);
1565                             }
1566
1567                             gmx_parallel_3dfft_execute(pfft_setup, GMX_FFT_REAL_TO_COMPLEX,
1568                                                        thread, wcycle);
1569                             if (thread == 0)
1570                             {
1571                                 wallcycle_stop(wcycle, ewcPME_FFT);
1572                             }
1573                         }
1574                     }
1575                     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1576                 }
1577                 bFirst = FALSE;
1578             }
1579             if (flags & GMX_PME_SOLVE)
1580             {
1581                 /* solve in k-space for our local cells */
1582 #pragma omp parallel num_threads(pme->nthread) private(thread)
1583                 {
1584                     try
1585                     {
1586                         int loop_count;
1587                         thread = gmx_omp_get_thread_num();
1588                         if (thread == 0)
1589                         {
1590                             wallcycle_start(wcycle, ewcLJPME);
1591                         }
1592
1593                         loop_count =
1594                             solve_pme_lj_yzx(pme, &pme->cfftgrid[2], TRUE,
1595                                              scaledBox[XX][XX]*scaledBox[YY][YY]*scaledBox[ZZ][ZZ],
1596                                              bCalcEnerVir,
1597                                              pme->nthread, thread);
1598                         if (thread == 0)
1599                         {
1600                             wallcycle_stop(wcycle, ewcLJPME);
1601                             inc_nrnb(nrnb, eNR_SOLVEPME, loop_count);
1602                         }
1603                     }
1604                     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1605                 }
1606             }
1607
1608             if (bCalcEnerVir)
1609             {
1610                 /* This should only be called on the master thread and
1611                  * after the threads have synchronized.
1612                  */
1613                 get_pme_ener_vir_lj(pme->solve_work, pme->nthread, &energy_AB[2+fep_state], vir_AB[2+fep_state]);
1614             }
1615
1616             if (bBackFFT)
1617             {
1618                 bFirst = !pme->doCoulomb;
1619                 calc_initial_lb_coeffs(pme, local_c6, local_sigma);
1620                 for (grid_index = 8; grid_index >= 2; --grid_index)
1621                 {
1622                     /* Unpack structure */
1623                     pmegrid    = &pme->pmegrid[grid_index];
1624                     fftgrid    = pme->fftgrid[grid_index];
1625                     pfft_setup = pme->pfft_setup[grid_index];
1626                     grid       = pmegrid->grid.grid;
1627                     calc_next_lb_coeffs(pme, local_sigma);
1628 #pragma omp parallel num_threads(pme->nthread) private(thread)
1629                     {
1630                         try
1631                         {
1632                             thread = gmx_omp_get_thread_num();
1633                             /* do 3d-invfft */
1634                             if (thread == 0)
1635                             {
1636                                 wallcycle_start(wcycle, ewcPME_FFT);
1637                             }
1638
1639                             gmx_parallel_3dfft_execute(pfft_setup, GMX_FFT_COMPLEX_TO_REAL,
1640                                                        thread, wcycle);
1641                             if (thread == 0)
1642                             {
1643                                 wallcycle_stop(wcycle, ewcPME_FFT);
1644
1645
1646                                 if (pme->nodeid == 0)
1647                                 {
1648                                     real ntot = pme->nkx*pme->nky*pme->nkz;
1649                                     npme  = static_cast<int>(ntot*std::log(ntot)/std::log(2.0));
1650                                     inc_nrnb(nrnb, eNR_FFT, 2*npme);
1651                                 }
1652                                 wallcycle_start(wcycle, ewcPME_GATHER);
1653                             }
1654
1655                             copy_fftgrid_to_pmegrid(pme, fftgrid, grid, grid_index, pme->nthread, thread);
1656                         }
1657                         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1658                     } /*#pragma omp parallel*/
1659
1660                     /* distribute local grid to all nodes */
1661 #if GMX_MPI
1662                     if (pme->nnodes > 1)
1663                     {
1664                         gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_BACKWARD);
1665                     }
1666 #endif
1667
1668                     unwrap_periodic_pmegrid(pme, grid);
1669
1670                     if (bCalcF)
1671                     {
1672                         /* interpolate forces for our local atoms */
1673                         bClearF = (bFirst && PAR(cr));
1674                         scale   = pme->bFEP ? (fep_state < 1 ? 1.0-lambda_lj : lambda_lj) : 1.0;
1675                         scale  *= lb_scale_factor[grid_index-2];
1676
1677 #pragma omp parallel for num_threads(pme->nthread) schedule(static)
1678                         for (thread = 0; thread < pme->nthread; thread++)
1679                         {
1680                             try
1681                             {
1682                                 gather_f_bsplines(pme, grid, bClearF, &pme->atc[0],
1683                                                   &pme->atc[0].spline[thread],
1684                                                   scale);
1685                             }
1686                             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1687                         }
1688
1689
1690                         inc_nrnb(nrnb, eNR_GATHERFBSP,
1691                                  pme->pme_order*pme->pme_order*pme->pme_order*pme->atc[0].n);
1692                     }
1693                     wallcycle_stop(wcycle, ewcPME_GATHER);
1694
1695                     bFirst = FALSE;
1696                 } /* for (grid_index = 8; grid_index >= 2; --grid_index) */
1697             }     /* if (bCalcF) */
1698         }         /* for (fep_state = 0; fep_state < fep_states_lj; ++fep_state) */
1699     }             /* if ((flags & GMX_PME_DO_LJ) && pme->ljpme_combination_rule == eljpmeLB) */
1700
1701     if (bCalcF && pme->nnodes > 1)
1702     {
1703         wallcycle_start(wcycle, ewcPME_REDISTXF);
1704         for (d = 0; d < pme->ndecompdim; d++)
1705         {
1706             atc = &pme->atc[d];
1707             if (d == pme->ndecompdim - 1)
1708             {
1709                 n_d = homenr;
1710                 f_d = f + start;
1711             }
1712             else
1713             {
1714                 n_d = pme->atc[d+1].n;
1715                 f_d = pme->atc[d+1].f;
1716             }
1717             if (DOMAINDECOMP(cr))
1718             {
1719                 dd_pmeredist_f(pme, atc, n_d, f_d,
1720                                d == pme->ndecompdim-1 && pme->bPPnode);
1721             }
1722         }
1723
1724         wallcycle_stop(wcycle, ewcPME_REDISTXF);
1725     }
1726
1727     if (bCalcEnerVir)
1728     {
1729         if (pme->doCoulomb)
1730         {
1731             if (!pme->bFEP_q)
1732             {
1733                 *energy_q = energy_AB[0];
1734                 m_add(vir_q, vir_AB[0], vir_q);
1735             }
1736             else
1737             {
1738                 *energy_q       = (1.0-lambda_q)*energy_AB[0] + lambda_q*energy_AB[1];
1739                 *dvdlambda_q   += energy_AB[1] - energy_AB[0];
1740                 for (i = 0; i < DIM; i++)
1741                 {
1742                     for (j = 0; j < DIM; j++)
1743                     {
1744                         vir_q[i][j] += (1.0-lambda_q)*vir_AB[0][i][j] +
1745                             lambda_q*vir_AB[1][i][j];
1746                     }
1747                 }
1748             }
1749             if (debug)
1750             {
1751                 fprintf(debug, "Electrostatic PME mesh energy: %g\n", *energy_q);
1752             }
1753         }
1754         else
1755         {
1756             *energy_q = 0;
1757         }
1758
1759         if (pme->doLJ)
1760         {
1761             if (!pme->bFEP_lj)
1762             {
1763                 *energy_lj = energy_AB[2];
1764                 m_add(vir_lj, vir_AB[2], vir_lj);
1765             }
1766             else
1767             {
1768                 *energy_lj     = (1.0-lambda_lj)*energy_AB[2] + lambda_lj*energy_AB[3];
1769                 *dvdlambda_lj += energy_AB[3] - energy_AB[2];
1770                 for (i = 0; i < DIM; i++)
1771                 {
1772                     for (j = 0; j < DIM; j++)
1773                     {
1774                         vir_lj[i][j] += (1.0-lambda_lj)*vir_AB[2][i][j] + lambda_lj*vir_AB[3][i][j];
1775                     }
1776                 }
1777             }
1778             if (debug)
1779             {
1780                 fprintf(debug, "Lennard-Jones PME mesh energy: %g\n", *energy_lj);
1781             }
1782         }
1783         else
1784         {
1785             *energy_lj = 0;
1786         }
1787     }
1788     return 0;
1789 }
1790
1791 void gmx_pme_destroy(gmx_pme_t *pme)
1792 {
1793     if (!pme)
1794     {
1795         return;
1796     }
1797
1798     delete pme->boxScaler;
1799
1800     sfree(pme->nnx);
1801     sfree(pme->nny);
1802     sfree(pme->nnz);
1803     sfree(pme->fshx);
1804     sfree(pme->fshy);
1805     sfree(pme->fshz);
1806
1807     for (int i = 0; i < pme->ngrids; ++i)
1808     {
1809         pmegrids_destroy(&pme->pmegrid[i]);
1810     }
1811     if (pme->pfft_setup)
1812     {
1813         for (int i = 0; i < pme->ngrids; ++i)
1814         {
1815             gmx_parallel_3dfft_destroy(pme->pfft_setup[i]);
1816         }
1817     }
1818     sfree(pme->fftgrid);
1819     sfree(pme->cfftgrid);
1820     sfree(pme->pfft_setup);
1821
1822     for (int i = 0; i < std::max(1, pme->ndecompdim); i++) //pme->atc[0] is always allocated
1823     {
1824         destroy_atomcomm(&pme->atc[i]);
1825     }
1826
1827     for (int i = 0; i < DIM; i++)
1828     {
1829         sfree(pme->bsp_mod[i]);
1830     }
1831
1832     sfree(pme->lb_buf1);
1833     sfree(pme->lb_buf2);
1834
1835     sfree(pme->bufv);
1836     sfree(pme->bufr);
1837
1838     if (pme->solve_work)
1839     {
1840         pme_free_all_work(&pme->solve_work, pme->nthread);
1841     }
1842
1843     sfree(pme->sum_qgrid_tmp);
1844     sfree(pme->sum_qgrid_dd_tmp);
1845
1846     destroy_pme_spline_work(pme->spline_work);
1847
1848     if (pme_gpu_active(pme) && pme->gpu)
1849     {
1850         pme_gpu_destroy(pme->gpu);
1851     }
1852
1853     delete pme;
1854 }
1855
1856 void gmx_pme_reinit_atoms(const gmx_pme_t *pme, const int nAtoms, const real *charges)
1857 {
1858     if (pme_gpu_active(pme))
1859     {
1860         pme_gpu_reinit_atoms(pme->gpu, nAtoms, charges);
1861     }
1862     // TODO: handle the CPU case here; handle the whole t_mdatoms
1863 }