src/gromacs/ewald/pme-only.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, 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 /* IMPORTANT FOR DEVELOPERS:
  38  *
  39  * Triclinic pme stuff isn't entirely trivial, and we've experienced
  40  * some bugs during development (many of them due to me). To avoid
  41  * this in the future, please check the following things if you make
  42  * changes in this file:
  43  *
  44  * 1. You should obtain identical (at least to the PME precision)
  45  *    energies, forces, and virial for
  46  *    a rectangular box and a triclinic one where the z (or y) axis is
  47  *    tilted a whole box side. For instance you could use these boxes:
  48  *
  49  *    rectangular       triclinic
  50  *     2  0  0           2  0  0
  51  *     0  2  0           0  2  0
  52  *     0  0  6           2  2  6
  53  *
  54  * 2. You should check the energy conservation in a triclinic box.
  55  *
  56  * It might seem an overkill, but better safe than sorry.
  57  * /Erik 001109
  58  */
  59
  60 #include "gmxpre.h"
  61
  62 #include <assert.h>
  63 #include <math.h>
  64 #include <stdio.h>
  65 #include <stdlib.h>
  66 #include <string.h>
  67
  68 #include "gromacs/ewald/pme.h"
  69 #include "gromacs/fft/parallel_3dfft.h"
  70 #include "gromacs/fileio/pdbio.h"
  71 #include "gromacs/gmxlib/network.h"
  72 #include "gromacs/gmxlib/nrnb.h"
  73 #include "gromacs/math/gmxcomplex.h"
  74 #include "gromacs/math/units.h"
  75 #include "gromacs/math/vec.h"
  76 #include "gromacs/mdtypes/commrec.h"
  77 #include "gromacs/mdtypes/inputrec.h"
  78 #include "gromacs/timing/cyclecounter.h"
  79 #include "gromacs/timing/wallcycle.h"
  80 #include "gromacs/utility/fatalerror.h"
  81 #include "gromacs/utility/futil.h"
  82 #include "gromacs/utility/gmxmpi.h"
  83 #include "gromacs/utility/gmxomp.h"
  84 #include "gromacs/utility/smalloc.h"
  85
  86 #include "pme-internal.h"
  87
  88 static void reset_pmeonly_counters(gmx_wallcycle_t wcycle,
  89                                    gmx_walltime_accounting_t walltime_accounting,
  90                                    t_nrnb *nrnb, t_inputrec *ir,
  91                                    gmx_int64_t step)
  92 {
  93     /* Reset all the counters related to performance over the run */
  94     wallcycle_stop(wcycle, ewcRUN);
  95     wallcycle_reset_all(wcycle);
  96     init_nrnb(nrnb);
  97     if (ir->nsteps >= 0)
  98     {
  99         /* ir->nsteps is not used here, but we update it for consistency */
 100         ir->nsteps -= step - ir->init_step;
 101     }
 102     ir->init_step = step;
 103     wallcycle_start(wcycle, ewcRUN);
 104     walltime_accounting_start(walltime_accounting);
 105 }
 106
 107
 108 static void gmx_pmeonly_switch(int *npmedata, struct gmx_pme_t ***pmedata,
 109                                ivec grid_size,
 110                                real ewaldcoeff_q, real ewaldcoeff_lj,
 111                                t_commrec *cr, t_inputrec *ir,
 112                                struct gmx_pme_t **pme_ret)
 113 {
 114     int               ind;
 115     struct gmx_pme_t *pme = nullptr;
 116
 117     ind = 0;
 118     while (ind < *npmedata)
 119     {
 120         pme = (*pmedata)[ind];
 121         if (pme->nkx == grid_size[XX] &&
 122             pme->nky == grid_size[YY] &&
 123             pme->nkz == grid_size[ZZ])
 124         {
 125             *pme_ret = pme;
 126
 127             return;
 128         }
 129
 130         ind++;
 131     }
 132
 133     (*npmedata)++;
 134     srenew(*pmedata, *npmedata);
 135
 136     /* Generate a new PME data structure, copying part of the old pointers */
 137     gmx_pme_reinit(&((*pmedata)[ind]), cr, pme, ir, grid_size, ewaldcoeff_q, ewaldcoeff_lj);
 138
 139     *pme_ret = (*pmedata)[ind];
 140 }
 141
 142 int gmx_pmeonly(struct gmx_pme_t *pme,
 143                 t_commrec *cr,    t_nrnb *mynrnb,
 144                 gmx_wallcycle_t wcycle,
 145                 gmx_walltime_accounting_t walltime_accounting,
 146                 real ewaldcoeff_q, real ewaldcoeff_lj,
 147                 t_inputrec *ir)
 148 {
 149     int                npmedata;
 150     struct gmx_pme_t **pmedata;
 151     gmx_pme_pp_t       pme_pp;
 152     int                ret;
 153     int                natoms;
 154     matrix             box;
 155     rvec              *x_pp       = nullptr, *f_pp = nullptr;
 156     real              *chargeA    = nullptr, *chargeB = nullptr;
 157     real              *c6A        = nullptr, *c6B = nullptr;
 158     real              *sigmaA     = nullptr, *sigmaB = nullptr;
 159     real               lambda_q   = 0;
 160     real               lambda_lj  = 0;
 161     int                maxshift_x = 0, maxshift_y = 0;
 162     real               energy_q, energy_lj, dvdlambda_q, dvdlambda_lj;
 163     matrix             vir_q, vir_lj;
 164     float              cycles;
 165     int                count;
 166     gmx_bool           bEnerVir;
 167     gmx_int64_t        step;
 168     ivec               grid_switch;
 169
 170     /* This data will only use with PME tuning, i.e. switching PME grids */
 171     npmedata = 1;
 172     snew(pmedata, npmedata);
 173     pmedata[0] = pme;
 174
 175     pme_pp = gmx_pme_pp_init(cr);
 176
 177     init_nrnb(mynrnb);
 178
 179     count = 0;
 180     do /****** this is a quasi-loop over time steps! */
 181     {
 182         /* The reason for having a loop here is PME grid tuning/switching */
 183         do
 184         {
 185             /* Domain decomposition */
 186             ret = gmx_pme_recv_coeffs_coords(pme_pp,
 187                                              &natoms,
 188                                              &chargeA, &chargeB,
 189                                              &c6A, &c6B,
 190                                              &sigmaA, &sigmaB,
 191                                              box, &x_pp, &f_pp,
 192                                              &maxshift_x, &maxshift_y,
 193                                              &lambda_q, &lambda_lj,
 194                                              &bEnerVir,
 195                                              &step,
 196                                              grid_switch, &ewaldcoeff_q, &ewaldcoeff_lj);
 197
 198             if (ret == pmerecvqxSWITCHGRID)
 199             {
 200                 /* Switch the PME grid to grid_switch */
 201                 gmx_pmeonly_switch(&npmedata, &pmedata, grid_switch, ewaldcoeff_q, ewaldcoeff_lj, cr, ir, &pme);
 202             }
 203
 204             if (ret == pmerecvqxRESETCOUNTERS)
 205             {
 206                 /* Reset the cycle and flop counters */
 207                 reset_pmeonly_counters(wcycle, walltime_accounting, mynrnb, ir, step);
 208             }
 209         }
 210         while (ret == pmerecvqxSWITCHGRID || ret == pmerecvqxRESETCOUNTERS);
 211
 212         if (ret == pmerecvqxFINISH)
 213         {
 214             /* We should stop: break out of the loop */
 215             break;
 216         }
 217
 218         if (count == 0)
 219         {
 220             wallcycle_start(wcycle, ewcRUN);
 221             walltime_accounting_start(walltime_accounting);
 222         }
 223
 224         wallcycle_start(wcycle, ewcPMEMESH);
 225
 226         dvdlambda_q  = 0;
 227         dvdlambda_lj = 0;
 228         clear_mat(vir_q);
 229         clear_mat(vir_lj);
 230
 231         gmx_pme_do(pme, 0, natoms, x_pp, f_pp,
 232                    chargeA, chargeB, c6A, c6B, sigmaA, sigmaB, box,
 233                    cr, maxshift_x, maxshift_y, mynrnb, wcycle,
 234                    vir_q, vir_lj,
 235                    &energy_q, &energy_lj, lambda_q, lambda_lj, &dvdlambda_q, &dvdlambda_lj,
 236                    GMX_PME_DO_ALL_F | (bEnerVir ? GMX_PME_CALC_ENER_VIR : 0));
 237
 238         cycles = wallcycle_stop(wcycle, ewcPMEMESH);
 239
 240         gmx_pme_send_force_vir_ener(pme_pp,
 241                                     f_pp, vir_q, energy_q, vir_lj, energy_lj,
 242                                     dvdlambda_q, dvdlambda_lj, cycles);
 243
 244         count++;
 245     } /***** end of quasi-loop, we stop with the break above */
 246     while (TRUE);
 247
 248     walltime_accounting_end(walltime_accounting);
 249
 250     return 0;
 251 }