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