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     real                 energy_AB[4];
1128     matrix               vir_AB[4];
1129     real                 scale, lambda;
1130     gmx_bool             bClearF;
1131     gmx_parallel_3dfft_t pfft_setup;
1132     real              *  fftgrid;
1133     t_complex          * cfftgrid;
1134     int                  thread;
1135     gmx_bool             bFirst, bDoSplines;
1136     int                  fep_state;
1137     int                  fep_states_lj           = pme->bFEP_lj ? 2 : 1;
1138     const gmx_bool       bCalcEnerVir            = (flags & GMX_PME_CALC_ENER_VIR) != 0;
1139     const gmx_bool       bBackFFT                = (flags & (GMX_PME_CALC_F | GMX_PME_CALC_POT)) != 0;
1140     const gmx_bool       bCalcF                  = (flags & GMX_PME_CALC_F) != 0;
1141
1142     /* We could be passing lambda!=1 while no q or LJ is actually perturbed */
1143     if (!pme->bFEP_q)
1144     {
1145         lambda_q  = 1;
1146     }
1147     if (!pme->bFEP_lj)
1148     {
1149         lambda_lj = 1;
1150     }
1151
1152     assert(pme->nnodes > 0);
1153     assert(pme->nnodes == 1 || pme->ndecompdim > 0);
1154
1155     if (pme->nnodes > 1)
1156     {
1157         atc      = &pme->atc[0];
1158         atc->npd = homenr;
1159         if (atc->npd > atc->pd_nalloc)
1160         {
1161             atc->pd_nalloc = over_alloc_dd(atc->npd);
1162             srenew(atc->pd, atc->pd_nalloc);
1163         }
1164         for (d = pme->ndecompdim-1; d >= 0; d--)
1165         {
1166             atc           = &pme->atc[d];
1167             atc->maxshift = (atc->dimind == 0 ? maxshift_x : maxshift_y);
1168         }
1169     }
1170     else
1171     {
1172         atc = &pme->atc[0];
1173         /* This could be necessary for TPI */
1174         pme->atc[0].n = homenr;
1175         if (DOMAINDECOMP(cr))
1176         {
1177             pme_realloc_atomcomm_things(atc);
1178         }
1179         atc->x = x;
1180         atc->f = f;
1181     }
1182
1183     matrix scaledBox;
1184     pme->boxScaler->scaleBox(box, scaledBox);
1185
1186     gmx::invertBoxMatrix(scaledBox, pme->recipbox);
1187     bFirst = TRUE;
1188
1189     /* For simplicity, we construct the splines for all particles if
1190      * more than one PME calculations is needed. Some optimization
1191      * could be done by keeping track of which atoms have splines
1192      * constructed, and construct new splines on each pass for atoms
1193      * that don't yet have them.
1194      */
1195
1196     bDoSplines = pme->bFEP || (pme->doCoulomb && pme->doLJ);
1197
1198     /* We need a maximum of four separate PME calculations:
1199      * grid_index=0: Coulomb PME with charges from state A
1200      * grid_index=1: Coulomb PME with charges from state B
1201      * grid_index=2: LJ PME with C6 from state A
1202      * grid_index=3: LJ PME with C6 from state B
1203      * For Lorentz-Berthelot combination rules, a separate loop is used to
1204      * calculate all the terms
1205      */
1206
1207     /* If we are doing LJ-PME with LB, we only do Q here */
1208     max_grid_index = (pme->ljpme_combination_rule == eljpmeLB) ? DO_Q : DO_Q_AND_LJ;
1209
1210     for (grid_index = 0; grid_index < max_grid_index; ++grid_index)
1211     {
1212         /* Check if we should do calculations at this grid_index
1213          * If grid_index is odd we should be doing FEP
1214          * If grid_index < 2 we should be doing electrostatic PME
1215          * If grid_index >= 2 we should be doing LJ-PME
1216          */
1217         if ((grid_index <  DO_Q && (!pme->doCoulomb ||
1218                                     (grid_index == 1 && !pme->bFEP_q))) ||
1219             (grid_index >= DO_Q && (!pme->doLJ ||
1220                                     (grid_index == 3 && !pme->bFEP_lj))))
1221         {
1222             continue;
1223         }
1224         /* Unpack structure */
1225         pmegrid    = &pme->pmegrid[grid_index];
1226         fftgrid    = pme->fftgrid[grid_index];
1227         cfftgrid   = pme->cfftgrid[grid_index];
1228         pfft_setup = pme->pfft_setup[grid_index];
1229         switch (grid_index)
1230         {
1231             case 0: coefficient = chargeA + start; break;
1232             case 1: coefficient = chargeB + start; break;
1233             case 2: coefficient = c6A + start; break;
1234             case 3: coefficient = c6B + start; break;
1235         }
1236
1237         grid = pmegrid->grid.grid;
1238
1239         if (debug)
1240         {
1241             fprintf(debug, "PME: number of ranks = %d, rank = %d\n",
1242                     cr->nnodes, cr->nodeid);
1243             fprintf(debug, "Grid = %p\n", static_cast<void*>(grid));
1244             if (grid == nullptr)
1245             {
1246                 gmx_fatal(FARGS, "No grid!");
1247             }
1248         }
1249
1250         if (pme->nnodes == 1)
1251         {
1252             atc->coefficient = coefficient;
1253         }
1254         else
1255         {
1256             wallcycle_start(wcycle, ewcPME_REDISTXF);
1257             do_redist_pos_coeffs(pme, cr, start, homenr, bFirst, x, coefficient);
1258
1259             wallcycle_stop(wcycle, ewcPME_REDISTXF);
1260         }
1261
1262         if (debug)
1263         {
1264             fprintf(debug, "Rank= %6d, pme local particles=%6d\n",
1265                     cr->nodeid, atc->n);
1266         }
1267
1268         if (flags & GMX_PME_SPREAD)
1269         {
1270             wallcycle_start(wcycle, ewcPME_SPREAD);
1271
1272             /* Spread the coefficients on a grid */
1273             spread_on_grid(pme, &pme->atc[0], pmegrid, bFirst, TRUE, fftgrid, bDoSplines, grid_index);
1274
1275             if (bFirst)
1276             {
1277                 inc_nrnb(nrnb, eNR_WEIGHTS, DIM*atc->n);
1278             }
1279             inc_nrnb(nrnb, eNR_SPREADBSP,
1280                      pme->pme_order*pme->pme_order*pme->pme_order*atc->n);
1281
1282             if (!pme->bUseThreads)
1283             {
1284                 wrap_periodic_pmegrid(pme, grid);
1285
1286                 /* sum contributions to local grid from other nodes */
1287 #if GMX_MPI
1288                 if (pme->nnodes > 1)
1289                 {
1290                     gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_FORWARD);
1291                 }
1292 #endif
1293
1294                 copy_pmegrid_to_fftgrid(pme, grid, fftgrid, grid_index);
1295             }
1296
1297             wallcycle_stop(wcycle, ewcPME_SPREAD);
1298
1299             /* TODO If the OpenMP and single-threaded implementations
1300                converge, then spread_on_grid() and
1301                copy_pmegrid_to_fftgrid() will perhaps live in the same
1302                source file.
1303              */
1304         }
1305
1306         /* Here we start a large thread parallel region */
1307 #pragma omp parallel num_threads(pme->nthread) private(thread)
1308         {
1309             try
1310             {
1311                 thread = gmx_omp_get_thread_num();
1312                 if (flags & GMX_PME_SOLVE)
1313                 {
1314                     int loop_count;
1315
1316                     /* do 3d-fft */
1317                     if (thread == 0)
1318                     {
1319                         wallcycle_start(wcycle, ewcPME_FFT);
1320                     }
1321                     gmx_parallel_3dfft_execute(pfft_setup, GMX_FFT_REAL_TO_COMPLEX,
1322                                                thread, wcycle);
1323                     if (thread == 0)
1324                     {
1325                         wallcycle_stop(wcycle, ewcPME_FFT);
1326                     }
1327
1328                     /* solve in k-space for our local cells */
1329                     if (thread == 0)
1330                     {
1331                         wallcycle_start(wcycle, (grid_index < DO_Q ? ewcPME_SOLVE : ewcLJPME));
1332                     }
1333                     if (grid_index < DO_Q)
1334                     {
1335                         loop_count =
1336                             solve_pme_yzx(pme, cfftgrid,
1337                                           scaledBox[XX][XX]*scaledBox[YY][YY]*scaledBox[ZZ][ZZ],
1338                                           bCalcEnerVir,
1339                                           pme->nthread, thread);
1340                     }
1341                     else
1342                     {
1343                         loop_count =
1344                             solve_pme_lj_yzx(pme, &cfftgrid, FALSE,
1345                                              scaledBox[XX][XX]*scaledBox[YY][YY]*scaledBox[ZZ][ZZ],
1346                                              bCalcEnerVir,
1347                                              pme->nthread, thread);
1348                     }
1349
1350                     if (thread == 0)
1351                     {
1352                         wallcycle_stop(wcycle, (grid_index < DO_Q ? ewcPME_SOLVE : ewcLJPME));
1353                         inc_nrnb(nrnb, eNR_SOLVEPME, loop_count);
1354                     }
1355                 }
1356
1357                 if (bBackFFT)
1358                 {
1359                     /* do 3d-invfft */
1360                     if (thread == 0)
1361                     {
1362                         wallcycle_start(wcycle, ewcPME_FFT);
1363                     }
1364                     gmx_parallel_3dfft_execute(pfft_setup, GMX_FFT_COMPLEX_TO_REAL,
1365                                                thread, wcycle);
1366                     if (thread == 0)
1367                     {
1368                         wallcycle_stop(wcycle, ewcPME_FFT);
1369
1370
1371                         if (pme->nodeid == 0)
1372                         {
1373                             real ntot = pme->nkx*pme->nky*pme->nkz;
1374                             npme  = static_cast<int>(ntot*std::log(ntot)/std::log(2.0));
1375                             inc_nrnb(nrnb, eNR_FFT, 2*npme);
1376                         }
1377
1378                         /* Note: this wallcycle region is closed below
1379                            outside an OpenMP region, so take care if
1380                            refactoring code here. */
1381                         wallcycle_start(wcycle, ewcPME_GATHER);
1382                     }
1383
1384                     copy_fftgrid_to_pmegrid(pme, fftgrid, grid, grid_index, pme->nthread, thread);
1385                 }
1386             } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1387         }
1388         /* End of thread parallel section.
1389          * With MPI we have to synchronize here before gmx_sum_qgrid_dd.
1390          */
1391
1392         if (bBackFFT)
1393         {
1394             /* distribute local grid to all nodes */
1395 #if GMX_MPI
1396             if (pme->nnodes > 1)
1397             {
1398                 gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_BACKWARD);
1399             }
1400 #endif
1401
1402             unwrap_periodic_pmegrid(pme, grid);
1403         }
1404
1405         if (bCalcF)
1406         {
1407             /* interpolate forces for our local atoms */
1408
1409
1410             /* If we are running without parallelization,
1411              * atc->f is the actual force array, not a buffer,
1412              * therefore we should not clear it.
1413              */
1414             lambda  = grid_index < DO_Q ? lambda_q : lambda_lj;
1415             bClearF = (bFirst && PAR(cr));
1416 #pragma omp parallel for num_threads(pme->nthread) schedule(static)
1417             for (thread = 0; thread < pme->nthread; thread++)
1418             {
1419                 try
1420                 {
1421                     gather_f_bsplines(pme, grid, bClearF, atc,
1422                                       &atc->spline[thread],
1423                                       pme->bFEP ? (grid_index % 2 == 0 ? 1.0-lambda : lambda) : 1.0);
1424                 }
1425                 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1426             }
1427
1428
1429             inc_nrnb(nrnb, eNR_GATHERFBSP,
1430                      pme->pme_order*pme->pme_order*pme->pme_order*pme->atc[0].n);
1431             /* Note: this wallcycle region is opened above inside an OpenMP
1432                region, so take care if refactoring code here. */
1433             wallcycle_stop(wcycle, ewcPME_GATHER);
1434         }
1435
1436         if (bCalcEnerVir)
1437         {
1438             /* This should only be called on the master thread
1439              * and after the threads have synchronized.
1440              */
1441             if (grid_index < 2)
1442             {
1443                 get_pme_ener_vir_q(pme->solve_work, pme->nthread, &energy_AB[grid_index], vir_AB[grid_index]);
1444             }
1445             else
1446             {
1447                 get_pme_ener_vir_lj(pme->solve_work, pme->nthread, &energy_AB[grid_index], vir_AB[grid_index]);
1448             }
1449         }
1450         bFirst = FALSE;
1451     } /* of grid_index-loop */
1452
1453     /* For Lorentz-Berthelot combination rules in LJ-PME, we need to calculate
1454      * seven terms. */
1455
1456     if (pme->doLJ && pme->ljpme_combination_rule == eljpmeLB)
1457     {
1458         /* Loop over A- and B-state if we are doing FEP */
1459         for (fep_state = 0; fep_state < fep_states_lj; ++fep_state)
1460         {
1461             real *local_c6 = nullptr, *local_sigma = nullptr, *RedistC6 = nullptr, *RedistSigma = nullptr;
1462             if (pme->nnodes == 1)
1463             {
1464                 if (pme->lb_buf1 == nullptr)
1465                 {
1466                     pme->lb_buf_nalloc = pme->atc[0].n;
1467                     snew(pme->lb_buf1, pme->lb_buf_nalloc);
1468                 }
1469                 pme->atc[0].coefficient = pme->lb_buf1;
1470                 switch (fep_state)
1471                 {
1472                     case 0:
1473                         local_c6      = c6A;
1474                         local_sigma   = sigmaA;
1475                         break;
1476                     case 1:
1477                         local_c6      = c6B;
1478                         local_sigma   = sigmaB;
1479                         break;
1480                     default:
1481                         gmx_incons("Trying to access wrong FEP-state in LJ-PME routine");
1482                 }
1483             }
1484             else
1485             {
1486                 atc = &pme->atc[0];
1487                 switch (fep_state)
1488                 {
1489                     case 0:
1490                         RedistC6      = c6A;
1491                         RedistSigma   = sigmaA;
1492                         break;
1493                     case 1:
1494                         RedistC6      = c6B;
1495                         RedistSigma   = sigmaB;
1496                         break;
1497                     default:
1498                         gmx_incons("Trying to access wrong FEP-state in LJ-PME routine");
1499                 }
1500                 wallcycle_start(wcycle, ewcPME_REDISTXF);
1501
1502                 do_redist_pos_coeffs(pme, cr, start, homenr, bFirst, x, RedistC6);
1503                 if (pme->lb_buf_nalloc < atc->n)
1504                 {
1505                     pme->lb_buf_nalloc = atc->nalloc;
1506                     srenew(pme->lb_buf1, pme->lb_buf_nalloc);
1507                     srenew(pme->lb_buf2, pme->lb_buf_nalloc);
1508                 }
1509                 local_c6 = pme->lb_buf1;
1510                 for (i = 0; i < atc->n; ++i)
1511                 {
1512                     local_c6[i] = atc->coefficient[i];
1513                 }
1514
1515                 do_redist_pos_coeffs(pme, cr, start, homenr, FALSE, x, RedistSigma);
1516                 local_sigma = pme->lb_buf2;
1517                 for (i = 0; i < atc->n; ++i)
1518                 {
1519                     local_sigma[i] = atc->coefficient[i];
1520                 }
1521
1522                 wallcycle_stop(wcycle, ewcPME_REDISTXF);
1523             }
1524             calc_initial_lb_coeffs(pme, local_c6, local_sigma);
1525
1526             /*Seven terms in LJ-PME with LB, grid_index < 2 reserved for electrostatics*/
1527             for (grid_index = 2; grid_index < 9; ++grid_index)
1528             {
1529                 /* Unpack structure */
1530                 pmegrid    = &pme->pmegrid[grid_index];
1531                 fftgrid    = pme->fftgrid[grid_index];
1532                 pfft_setup = pme->pfft_setup[grid_index];
1533                 calc_next_lb_coeffs(pme, local_sigma);
1534                 grid = pmegrid->grid.grid;
1535
1536                 if (flags & GMX_PME_SPREAD)
1537                 {
1538                     wallcycle_start(wcycle, ewcPME_SPREAD);
1539                     /* Spread the c6 on a grid */
1540                     spread_on_grid(pme, &pme->atc[0], pmegrid, bFirst, TRUE, fftgrid, bDoSplines, grid_index);
1541
1542                     if (bFirst)
1543                     {
1544                         inc_nrnb(nrnb, eNR_WEIGHTS, DIM*atc->n);
1545                     }
1546
1547                     inc_nrnb(nrnb, eNR_SPREADBSP,
1548                              pme->pme_order*pme->pme_order*pme->pme_order*atc->n);
1549                     if (pme->nthread == 1)
1550                     {
1551                         wrap_periodic_pmegrid(pme, grid);
1552                         /* sum contributions to local grid from other nodes */
1553 #if GMX_MPI
1554                         if (pme->nnodes > 1)
1555                         {
1556                             gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_FORWARD);
1557                         }
1558 #endif
1559                         copy_pmegrid_to_fftgrid(pme, grid, fftgrid, grid_index);
1560                     }
1561                     wallcycle_stop(wcycle, ewcPME_SPREAD);
1562                 }
1563                 /*Here we start a large thread parallel region*/
1564 #pragma omp parallel num_threads(pme->nthread) private(thread)
1565                 {
1566                     try
1567                     {
1568                         thread = gmx_omp_get_thread_num();
1569                         if (flags & GMX_PME_SOLVE)
1570                         {
1571                             /* do 3d-fft */
1572                             if (thread == 0)
1573                             {
1574                                 wallcycle_start(wcycle, ewcPME_FFT);
1575                             }
1576
1577                             gmx_parallel_3dfft_execute(pfft_setup, GMX_FFT_REAL_TO_COMPLEX,
1578                                                        thread, wcycle);
1579                             if (thread == 0)
1580                             {
1581                                 wallcycle_stop(wcycle, ewcPME_FFT);
1582                             }
1583                         }
1584                     }
1585                     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1586                 }
1587                 bFirst = FALSE;
1588             }
1589             if (flags & GMX_PME_SOLVE)
1590             {
1591                 /* solve in k-space for our local cells */
1592 #pragma omp parallel num_threads(pme->nthread) private(thread)
1593                 {
1594                     try
1595                     {
1596                         int loop_count;
1597                         thread = gmx_omp_get_thread_num();
1598                         if (thread == 0)
1599                         {
1600                             wallcycle_start(wcycle, ewcLJPME);
1601                         }
1602
1603                         loop_count =
1604                             solve_pme_lj_yzx(pme, &pme->cfftgrid[2], TRUE,
1605                                              scaledBox[XX][XX]*scaledBox[YY][YY]*scaledBox[ZZ][ZZ],
1606                                              bCalcEnerVir,
1607                                              pme->nthread, thread);
1608                         if (thread == 0)
1609                         {
1610                             wallcycle_stop(wcycle, ewcLJPME);
1611                             inc_nrnb(nrnb, eNR_SOLVEPME, loop_count);
1612                         }
1613                     }
1614                     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1615                 }
1616             }
1617
1618             if (bCalcEnerVir)
1619             {
1620                 /* This should only be called on the master thread and
1621                  * after the threads have synchronized.
1622                  */
1623                 get_pme_ener_vir_lj(pme->solve_work, pme->nthread, &energy_AB[2+fep_state], vir_AB[2+fep_state]);
1624             }
1625
1626             if (bBackFFT)
1627             {
1628                 bFirst = !pme->doCoulomb;
1629                 calc_initial_lb_coeffs(pme, local_c6, local_sigma);
1630                 for (grid_index = 8; grid_index >= 2; --grid_index)
1631                 {
1632                     /* Unpack structure */
1633                     pmegrid    = &pme->pmegrid[grid_index];
1634                     fftgrid    = pme->fftgrid[grid_index];
1635                     pfft_setup = pme->pfft_setup[grid_index];
1636                     grid       = pmegrid->grid.grid;
1637                     calc_next_lb_coeffs(pme, local_sigma);
1638 #pragma omp parallel num_threads(pme->nthread) private(thread)
1639                     {
1640                         try
1641                         {
1642                             thread = gmx_omp_get_thread_num();
1643                             /* do 3d-invfft */
1644                             if (thread == 0)
1645                             {
1646                                 wallcycle_start(wcycle, ewcPME_FFT);
1647                             }
1648
1649                             gmx_parallel_3dfft_execute(pfft_setup, GMX_FFT_COMPLEX_TO_REAL,
1650                                                        thread, wcycle);
1651                             if (thread == 0)
1652                             {
1653                                 wallcycle_stop(wcycle, ewcPME_FFT);
1654
1655
1656                                 if (pme->nodeid == 0)
1657                                 {
1658                                     real ntot = pme->nkx*pme->nky*pme->nkz;
1659                                     npme  = static_cast<int>(ntot*std::log(ntot)/std::log(2.0));
1660                                     inc_nrnb(nrnb, eNR_FFT, 2*npme);
1661                                 }
1662                                 wallcycle_start(wcycle, ewcPME_GATHER);
1663                             }
1664
1665                             copy_fftgrid_to_pmegrid(pme, fftgrid, grid, grid_index, pme->nthread, thread);
1666                         }
1667                         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1668                     } /*#pragma omp parallel*/
1669
1670                     /* distribute local grid to all nodes */
1671 #if GMX_MPI
1672                     if (pme->nnodes > 1)
1673                     {
1674                         gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_BACKWARD);
1675                     }
1676 #endif
1677
1678                     unwrap_periodic_pmegrid(pme, grid);
1679
1680                     if (bCalcF)
1681                     {
1682                         /* interpolate forces for our local atoms */
1683                         bClearF = (bFirst && PAR(cr));
1684                         scale   = pme->bFEP ? (fep_state < 1 ? 1.0-lambda_lj : lambda_lj) : 1.0;
1685                         scale  *= lb_scale_factor[grid_index-2];
1686
1687 #pragma omp parallel for num_threads(pme->nthread) schedule(static)
1688                         for (thread = 0; thread < pme->nthread; thread++)
1689                         {
1690                             try
1691                             {
1692                                 gather_f_bsplines(pme, grid, bClearF, &pme->atc[0],
1693                                                   &pme->atc[0].spline[thread],
1694                                                   scale);
1695                             }
1696                             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1697                         }
1698
1699
1700                         inc_nrnb(nrnb, eNR_GATHERFBSP,
1701                                  pme->pme_order*pme->pme_order*pme->pme_order*pme->atc[0].n);
1702                     }
1703                     wallcycle_stop(wcycle, ewcPME_GATHER);
1704
1705                     bFirst = FALSE;
1706                 } /* for (grid_index = 8; grid_index >= 2; --grid_index) */
1707             }     /* if (bCalcF) */
1708         }         /* for (fep_state = 0; fep_state < fep_states_lj; ++fep_state) */
1709     }             /* if ((flags & GMX_PME_DO_LJ) && pme->ljpme_combination_rule == eljpmeLB) */
1710
1711     if (bCalcF && pme->nnodes > 1)
1712     {
1713         wallcycle_start(wcycle, ewcPME_REDISTXF);
1714         for (d = 0; d < pme->ndecompdim; d++)
1715         {
1716             atc = &pme->atc[d];
1717             if (d == pme->ndecompdim - 1)
1718             {
1719                 n_d = homenr;
1720                 f_d = f + start;
1721             }
1722             else
1723             {
1724                 n_d = pme->atc[d+1].n;
1725                 f_d = pme->atc[d+1].f;
1726             }
1727             if (DOMAINDECOMP(cr))
1728             {
1729                 dd_pmeredist_f(pme, atc, n_d, f_d,
1730                                d == pme->ndecompdim-1 && pme->bPPnode);
1731             }
1732         }
1733
1734         wallcycle_stop(wcycle, ewcPME_REDISTXF);
1735     }
1736
1737     if (bCalcEnerVir)
1738     {
1739         if (pme->doCoulomb)
1740         {
1741             if (!pme->bFEP_q)
1742             {
1743                 *energy_q = energy_AB[0];
1744                 m_add(vir_q, vir_AB[0], vir_q);
1745             }
1746             else
1747             {
1748                 *energy_q       = (1.0-lambda_q)*energy_AB[0] + lambda_q*energy_AB[1];
1749                 *dvdlambda_q   += energy_AB[1] - energy_AB[0];
1750                 for (i = 0; i < DIM; i++)
1751                 {
1752                     for (j = 0; j < DIM; j++)
1753                     {
1754                         vir_q[i][j] += (1.0-lambda_q)*vir_AB[0][i][j] +
1755                             lambda_q*vir_AB[1][i][j];
1756                     }
1757                 }
1758             }
1759             if (debug)
1760             {
1761                 fprintf(debug, "Electrostatic PME mesh energy: %g\n", *energy_q);
1762             }
1763         }
1764         else
1765         {
1766             *energy_q = 0;
1767         }
1768
1769         if (pme->doLJ)
1770         {
1771             if (!pme->bFEP_lj)
1772             {
1773                 *energy_lj = energy_AB[2];
1774                 m_add(vir_lj, vir_AB[2], vir_lj);
1775             }
1776             else
1777             {
1778                 *energy_lj     = (1.0-lambda_lj)*energy_AB[2] + lambda_lj*energy_AB[3];
1779                 *dvdlambda_lj += energy_AB[3] - energy_AB[2];
1780                 for (i = 0; i < DIM; i++)
1781                 {
1782                     for (j = 0; j < DIM; j++)
1783                     {
1784                         vir_lj[i][j] += (1.0-lambda_lj)*vir_AB[2][i][j] + lambda_lj*vir_AB[3][i][j];
1785                     }
1786                 }
1787             }
1788             if (debug)
1789             {
1790                 fprintf(debug, "Lennard-Jones PME mesh energy: %g\n", *energy_lj);
1791             }
1792         }
1793         else
1794         {
1795             *energy_lj = 0;
1796         }
1797     }
1798     return 0;
1799 }
1800
1801 void gmx_pme_destroy(gmx_pme_t *pme)
1802 {
1803     if (!pme)
1804     {
1805         return;
1806     }
1807
1808     delete pme->boxScaler;
1809
1810     sfree(pme->nnx);
1811     sfree(pme->nny);
1812     sfree(pme->nnz);
1813     sfree(pme->fshx);
1814     sfree(pme->fshy);
1815     sfree(pme->fshz);
1816
1817     for (int i = 0; i < pme->ngrids; ++i)
1818     {
1819         pmegrids_destroy(&pme->pmegrid[i]);
1820     }
1821     if (pme->pfft_setup)
1822     {
1823         for (int i = 0; i < pme->ngrids; ++i)
1824         {
1825             gmx_parallel_3dfft_destroy(pme->pfft_setup[i]);
1826         }
1827     }
1828     sfree(pme->fftgrid);
1829     sfree(pme->cfftgrid);
1830     sfree(pme->pfft_setup);
1831
1832     for (int i = 0; i < std::max(1, pme->ndecompdim); i++) //pme->atc[0] is always allocated
1833     {
1834         destroy_atomcomm(&pme->atc[i]);
1835     }
1836
1837     for (int i = 0; i < DIM; i++)
1838     {
1839         sfree(pme->bsp_mod[i]);
1840     }
1841
1842     sfree(pme->lb_buf1);
1843     sfree(pme->lb_buf2);
1844
1845     sfree(pme->bufv);
1846     sfree(pme->bufr);
1847
1848     if (pme->solve_work)
1849     {
1850         pme_free_all_work(&pme->solve_work, pme->nthread);
1851     }
1852
1853     sfree(pme->sum_qgrid_tmp);
1854     sfree(pme->sum_qgrid_dd_tmp);
1855
1856     destroy_pme_spline_work(pme->spline_work);
1857
1858     if (pme_gpu_active(pme) && pme->gpu)
1859     {
1860         pme_gpu_destroy(pme->gpu);
1861     }
1862
1863     delete pme;
1864 }
1865
1866 void gmx_pme_reinit_atoms(const gmx_pme_t *pme, const int nAtoms, const real *charges)
1867 {
1868     if (pme_gpu_active(pme))
1869     {
1870         pme_gpu_reinit_atoms(pme->gpu, nAtoms, charges);
1871     }
1872     // TODO: handle the CPU case here; handle the whole t_mdatoms
1873 }