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, 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                                t_commrec *cr, t_inputrec *ir,
 111                                struct gmx_pme_t **pme_ret)
 112 {
 113     int               ind;
 114     struct gmx_pme_t *pme = NULL;
 115
 116     ind = 0;
 117     while (ind < *npmedata)
 118     {
 119         pme = (*pmedata)[ind];
 120         if (pme->nkx == grid_size[XX] &&
 121             pme->nky == grid_size[YY] &&
 122             pme->nkz == grid_size[ZZ])
 123         {
 124             *pme_ret = pme;
 125
 126             return;
 127         }
 128
 129         ind++;
 130     }
 131
 132     (*npmedata)++;
 133     srenew(*pmedata, *npmedata);
 134
 135     /* Generate a new PME data structure, copying part of the old pointers */
 136     gmx_pme_reinit(&((*pmedata)[ind]), cr, pme, ir, grid_size);
 137
 138     *pme_ret = (*pmedata)[ind];
 139 }
 140
 141 int gmx_pmeonly(struct gmx_pme_t *pme,
 142                 t_commrec *cr,    t_nrnb *mynrnb,
 143                 gmx_wallcycle_t wcycle,
 144                 gmx_walltime_accounting_t walltime_accounting,
 145                 real ewaldcoeff_q, real ewaldcoeff_lj,
 146                 t_inputrec *ir)
 147 {
 148     int                npmedata;
 149     struct gmx_pme_t **pmedata;
 150     gmx_pme_pp_t       pme_pp;
 151     int                ret;
 152     int                natoms;
 153     matrix             box;
 154     rvec              *x_pp       = NULL, *f_pp = NULL;
 155     real              *chargeA    = NULL, *chargeB = NULL;
 156     real              *c6A        = NULL, *c6B = NULL;
 157     real              *sigmaA     = NULL, *sigmaB = NULL;
 158     real               lambda_q   = 0;
 159     real               lambda_lj  = 0;
 160     int                maxshift_x = 0, maxshift_y = 0;
 161     real               energy_q, energy_lj, dvdlambda_q, dvdlambda_lj;
 162     matrix             vir_q, vir_lj;
 163     float              cycles;
 164     int                count;
 165     gmx_bool           bEnerVir;
 166     gmx_int64_t        step;
 167     ivec               grid_switch;
 168
 169     /* This data will only use with PME tuning, i.e. switching PME grids */
 170     npmedata = 1;
 171     snew(pmedata, npmedata);
 172     pmedata[0] = pme;
 173
 174     pme_pp = gmx_pme_pp_init(cr);
 175
 176     init_nrnb(mynrnb);
 177
 178     count = 0;
 179     do /****** this is a quasi-loop over time steps! */
 180     {
 181         /* The reason for having a loop here is PME grid tuning/switching */
 182         do
 183         {
 184             /* Domain decomposition */
 185             ret = gmx_pme_recv_coeffs_coords(pme_pp,
 186                                              &natoms,
 187                                              &chargeA, &chargeB,
 188                                              &c6A, &c6B,
 189                                              &sigmaA, &sigmaB,
 190                                              box, &x_pp, &f_pp,
 191                                              &maxshift_x, &maxshift_y,
 192                                              &lambda_q, &lambda_lj,
 193                                              &bEnerVir,
 194                                              &step,
 195                                              grid_switch, &ewaldcoeff_q, &ewaldcoeff_lj);
 196
 197             if (ret == pmerecvqxSWITCHGRID)
 198             {
 199                 /* Switch the PME grid to grid_switch */
 200                 gmx_pmeonly_switch(&npmedata, &pmedata, grid_switch, cr, ir, &pme);
 201             }
 202
 203             if (ret == pmerecvqxRESETCOUNTERS)
 204             {
 205                 /* Reset the cycle and flop counters */
 206                 reset_pmeonly_counters(wcycle, walltime_accounting, mynrnb, ir, step);
 207             }
 208         }
 209         while (ret == pmerecvqxSWITCHGRID || ret == pmerecvqxRESETCOUNTERS);
 210
 211         if (ret == pmerecvqxFINISH)
 212         {
 213             /* We should stop: break out of the loop */
 214             break;
 215         }
 216
 217         if (count == 0)
 218         {
 219             wallcycle_start(wcycle, ewcRUN);
 220             walltime_accounting_start(walltime_accounting);
 221         }
 222
 223         wallcycle_start(wcycle, ewcPMEMESH);
 224
 225         dvdlambda_q  = 0;
 226         dvdlambda_lj = 0;
 227         clear_mat(vir_q);
 228         clear_mat(vir_lj);
 229
 230         gmx_pme_do(pme, 0, natoms, x_pp, f_pp,
 231                    chargeA, chargeB, c6A, c6B, sigmaA, sigmaB, box,
 232                    cr, maxshift_x, maxshift_y, mynrnb, wcycle,
 233                    vir_q, ewaldcoeff_q, vir_lj, ewaldcoeff_lj,
 234                    &energy_q, &energy_lj, lambda_q, lambda_lj, &dvdlambda_q, &dvdlambda_lj,
 235                    GMX_PME_DO_ALL_F | (bEnerVir ? GMX_PME_CALC_ENER_VIR : 0));
 236
 237         cycles = wallcycle_stop(wcycle, ewcPMEMESH);
 238
 239         gmx_pme_send_force_vir_ener(pme_pp,
 240                                     f_pp, vir_q, energy_q, vir_lj, energy_lj,
 241                                     dvdlambda_q, dvdlambda_lj, cycles);
 242
 243         count++;
 244     } /***** end of quasi-loop, we stop with the break above */
 245     while (TRUE);
 246
 247     walltime_accounting_end(walltime_accounting);
 248
 249     return 0;
 250 }