src/gromacs/ewald/pme.cpp

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