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