src/gromacs/mdlib/stat.c

   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, 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 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <string.h>
  42 #include <stdio.h>
  43 #include "typedefs.h"
  44 #include "sysstuff.h"
  45 #include "gmx_fatal.h"
  46 #include "network.h"
  47 #include "txtdump.h"
  48 #include "names.h"
  49 #include "physics.h"
  50 #include "vec.h"
  51 #include "gromacs/math/utilities.h"
  52 #include "main.h"
  53 #include "force.h"
  54 #include "vcm.h"
  55 #include "smalloc.h"
  56 #include "gromacs/fileio/futil.h"
  57 #include "network.h"
  58 #include "rbin.h"
  59 #include "tgroup.h"
  60 #include "gromacs/fileio/xtcio.h"
  61 #include "gromacs/fileio/gmxfio.h"
  62 #include "gromacs/fileio/trnio.h"
  63 #include "domdec.h"
  64 #include "constr.h"
  65 #include "checkpoint.h"
  66 #include "xvgr.h"
  67 #include "md_support.h"
  68 #include "mdrun.h"
  69 #include "sim_util.h"
  70
  71 typedef struct gmx_global_stat
  72 {
  73     t_bin *rb;
  74     int   *itc0;
  75     int   *itc1;
  76 } t_gmx_global_stat;
  77
  78 gmx_global_stat_t global_stat_init(t_inputrec *ir)
  79 {
  80     gmx_global_stat_t gs;
  81
  82     snew(gs, 1);
  83
  84     gs->rb = mk_bin();
  85     snew(gs->itc0, ir->opts.ngtc);
  86     snew(gs->itc1, ir->opts.ngtc);
  87
  88     return gs;
  89 }
  90
  91 void global_stat_destroy(gmx_global_stat_t gs)
  92 {
  93     destroy_bin(gs->rb);
  94     sfree(gs->itc0);
  95     sfree(gs->itc1);
  96     sfree(gs);
  97 }
  98
  99 static int filter_enerdterm(real *afrom, gmx_bool bToBuffer, real *ato,
 100                             gmx_bool bTemp, gmx_bool bPres, gmx_bool bEner)
 101 {
 102     int i, to, from;
 103
 104     from = 0;
 105     to   = 0;
 106     for (i = 0; i < F_NRE; i++)
 107     {
 108         if (bToBuffer)
 109         {
 110             from = i;
 111         }
 112         else
 113         {
 114             to = i;
 115         }
 116         switch (i)
 117         {
 118             case F_EKIN:
 119             case F_TEMP:
 120             case F_DKDL:
 121                 if (bTemp)
 122                 {
 123                     ato[to++] = afrom[from++];
 124                 }
 125                 break;
 126             case F_PRES:
 127             case F_PDISPCORR:
 128                 if (bPres)
 129                 {
 130                     ato[to++] = afrom[from++];
 131                 }
 132                 break;
 133             default:
 134                 if (bEner)
 135                 {
 136                     ato[to++] = afrom[from++];
 137                 }
 138                 break;
 139         }
 140     }
 141
 142     return to;
 143 }
 144
 145 void global_stat(FILE *fplog, gmx_global_stat_t gs,
 146                  t_commrec *cr, gmx_enerdata_t *enerd,
 147                  tensor fvir, tensor svir, rvec mu_tot,
 148                  t_inputrec *inputrec,
 149                  gmx_ekindata_t *ekind, gmx_constr_t constr,
 150                  t_vcm *vcm,
 151                  int nsig, real *sig,
 152                  gmx_mtop_t *top_global, t_state *state_local,
 153                  gmx_bool bSumEkinhOld, int flags)
 154 /* instead of current system, gmx_booleans for summing virial, kinetic energy, and other terms */
 155 {
 156     t_bin     *rb;
 157     int       *itc0, *itc1;
 158     int        ie    = 0, ifv = 0, isv = 0, irmsd = 0, imu = 0;
 159     int        idedl = 0, idvdll = 0, idvdlnl = 0, iepl = 0, icm = 0, imass = 0, ica = 0, inb = 0;
 160     int        isig  = -1;
 161     int        icj   = -1, ici = -1, icx = -1;
 162     int        inn[egNR];
 163     real       copyenerd[F_NRE];
 164     int        nener, j;
 165     real      *rmsd_data = NULL;
 166     double     nb;
 167     gmx_bool   bVV, bTemp, bEner, bPres, bConstrVir, bEkinAveVel, bFirstIterate, bReadEkin;
 168
 169     bVV           = EI_VV(inputrec->eI);
 170     bTemp         = flags & CGLO_TEMPERATURE;
 171     bEner         = flags & CGLO_ENERGY;
 172     bPres         = (flags & CGLO_PRESSURE);
 173     bConstrVir    = (flags & CGLO_CONSTRAINT);
 174     bFirstIterate = (flags & CGLO_FIRSTITERATE);
 175     bEkinAveVel   = (inputrec->eI == eiVV || (inputrec->eI == eiVVAK && bPres));
 176     bReadEkin     = (flags & CGLO_READEKIN);
 177
 178     rb   = gs->rb;
 179     itc0 = gs->itc0;
 180     itc1 = gs->itc1;
 181
 182
 183     reset_bin(rb);
 184     /* This routine copies all the data to be summed to one big buffer
 185      * using the t_bin struct.
 186      */
 187
 188     /* First, we neeed to identify which enerd->term should be
 189        communicated.  Temperature and pressure terms should only be
 190        communicated and summed when they need to be, to avoid repeating
 191        the sums and overcounting. */
 192
 193     nener = filter_enerdterm(enerd->term, TRUE, copyenerd, bTemp, bPres, bEner);
 194
 195     /* First, the data that needs to be communicated with velocity verlet every time
 196        This is just the constraint virial.*/
 197     if (bConstrVir)
 198     {
 199         isv = add_binr(rb, DIM*DIM, svir[0]);
 200         where();
 201     }
 202
 203 /* We need the force virial and the kinetic energy for the first time through with velocity verlet */
 204     if (bTemp || !bVV)
 205     {
 206         if (ekind)
 207         {
 208             for (j = 0; (j < inputrec->opts.ngtc); j++)
 209             {
 210                 if (bSumEkinhOld)
 211                 {
 212                     itc0[j] = add_binr(rb, DIM*DIM, ekind->tcstat[j].ekinh_old[0]);
 213                 }
 214                 if (bEkinAveVel && !bReadEkin)
 215                 {
 216                     itc1[j] = add_binr(rb, DIM*DIM, ekind->tcstat[j].ekinf[0]);
 217                 }
 218                 else if (!bReadEkin)
 219                 {
 220                     itc1[j] = add_binr(rb, DIM*DIM, ekind->tcstat[j].ekinh[0]);
 221                 }
 222             }
 223             /* these probably need to be put into one of these categories */
 224             where();
 225             idedl = add_binr(rb, 1, &(ekind->dekindl));
 226             where();
 227             ica   = add_binr(rb, 1, &(ekind->cosacc.mvcos));
 228             where();
 229         }
 230     }
 231     where();
 232
 233     if ((bPres || !bVV) && bFirstIterate)
 234     {
 235         ifv = add_binr(rb, DIM*DIM, fvir[0]);
 236     }
 237
 238
 239     if (bEner)
 240     {
 241         where();
 242         if (bFirstIterate)
 243         {
 244             ie  = add_binr(rb, nener, copyenerd);
 245         }
 246         where();
 247         if (constr)
 248         {
 249             rmsd_data = constr_rmsd_data(constr);
 250             if (rmsd_data)
 251             {
 252                 irmsd = add_binr(rb, inputrec->eI == eiSD2 ? 3 : 2, rmsd_data);
 253             }
 254         }
 255         if (!NEED_MUTOT(*inputrec))
 256         {
 257             imu = add_binr(rb, DIM, mu_tot);
 258             where();
 259         }
 260
 261         if (bFirstIterate)
 262         {
 263             for (j = 0; (j < egNR); j++)
 264             {
 265                 inn[j] = add_binr(rb, enerd->grpp.nener, enerd->grpp.ener[j]);
 266             }
 267             where();
 268             if (inputrec->efep != efepNO)
 269             {
 270                 idvdll  = add_bind(rb, efptNR, enerd->dvdl_lin);
 271                 idvdlnl = add_bind(rb, efptNR, enerd->dvdl_nonlin);
 272                 if (enerd->n_lambda > 0)
 273                 {
 274                     iepl = add_bind(rb, enerd->n_lambda, enerd->enerpart_lambda);
 275                 }
 276             }
 277         }
 278     }
 279
 280     if (vcm)
 281     {
 282         icm   = add_binr(rb, DIM*vcm->nr, vcm->group_p[0]);
 283         where();
 284         imass = add_binr(rb, vcm->nr, vcm->group_mass);
 285         where();
 286         if (vcm->mode == ecmANGULAR)
 287         {
 288             icj   = add_binr(rb, DIM*vcm->nr, vcm->group_j[0]);
 289             where();
 290             icx   = add_binr(rb, DIM*vcm->nr, vcm->group_x[0]);
 291             where();
 292             ici   = add_binr(rb, DIM*DIM*vcm->nr, vcm->group_i[0][0]);
 293             where();
 294         }
 295     }
 296
 297     if (DOMAINDECOMP(cr))
 298     {
 299         nb  = cr->dd->nbonded_local;
 300         inb = add_bind(rb, 1, &nb);
 301     }
 302     where();
 303     if (nsig > 0)
 304     {
 305         isig = add_binr(rb, nsig, sig);
 306     }
 307
 308     /* Global sum it all */
 309     if (debug)
 310     {
 311         fprintf(debug, "Summing %d energies\n", rb->maxreal);
 312     }
 313     sum_bin(rb, cr);
 314     where();
 315
 316     /* Extract all the data locally */
 317
 318     if (bConstrVir)
 319     {
 320         extract_binr(rb, isv, DIM*DIM, svir[0]);
 321     }
 322
 323     /* We need the force virial and the kinetic energy for the first time through with velocity verlet */
 324     if (bTemp || !bVV)
 325     {
 326         if (ekind)
 327         {
 328             for (j = 0; (j < inputrec->opts.ngtc); j++)
 329             {
 330                 if (bSumEkinhOld)
 331                 {
 332                     extract_binr(rb, itc0[j], DIM*DIM, ekind->tcstat[j].ekinh_old[0]);
 333                 }
 334                 if (bEkinAveVel && !bReadEkin)
 335                 {
 336                     extract_binr(rb, itc1[j], DIM*DIM, ekind->tcstat[j].ekinf[0]);
 337                 }
 338                 else if (!bReadEkin)
 339                 {
 340                     extract_binr(rb, itc1[j], DIM*DIM, ekind->tcstat[j].ekinh[0]);
 341                 }
 342             }
 343             extract_binr(rb, idedl, 1, &(ekind->dekindl));
 344             extract_binr(rb, ica, 1, &(ekind->cosacc.mvcos));
 345             where();
 346         }
 347     }
 348     if ((bPres || !bVV) && bFirstIterate)
 349     {
 350         extract_binr(rb, ifv, DIM*DIM, fvir[0]);
 351     }
 352
 353     if (bEner)
 354     {
 355         if (bFirstIterate)
 356         {
 357             extract_binr(rb, ie, nener, copyenerd);
 358             if (rmsd_data)
 359             {
 360                 extract_binr(rb, irmsd, inputrec->eI == eiSD2 ? 3 : 2, rmsd_data);
 361             }
 362             if (!NEED_MUTOT(*inputrec))
 363             {
 364                 extract_binr(rb, imu, DIM, mu_tot);
 365             }
 366
 367             for (j = 0; (j < egNR); j++)
 368             {
 369                 extract_binr(rb, inn[j], enerd->grpp.nener, enerd->grpp.ener[j]);
 370             }
 371             if (inputrec->efep != efepNO)
 372             {
 373                 extract_bind(rb, idvdll, efptNR, enerd->dvdl_lin);
 374                 extract_bind(rb, idvdlnl, efptNR, enerd->dvdl_nonlin);
 375                 if (enerd->n_lambda > 0)
 376                 {
 377                     extract_bind(rb, iepl, enerd->n_lambda, enerd->enerpart_lambda);
 378                 }
 379             }
 380             if (DOMAINDECOMP(cr))
 381             {
 382                 extract_bind(rb, inb, 1, &nb);
 383                 if ((int)(nb + 0.5) != cr->dd->nbonded_global)
 384                 {
 385                     dd_print_missing_interactions(fplog, cr, (int)(nb + 0.5), top_global, state_local);
 386                 }
 387             }
 388             where();
 389
 390             filter_enerdterm(copyenerd, FALSE, enerd->term, bTemp, bPres, bEner);
 391         }
 392     }
 393
 394     if (vcm)
 395     {
 396         extract_binr(rb, icm, DIM*vcm->nr, vcm->group_p[0]);
 397         where();
 398         extract_binr(rb, imass, vcm->nr, vcm->group_mass);
 399         where();
 400         if (vcm->mode == ecmANGULAR)
 401         {
 402             extract_binr(rb, icj, DIM*vcm->nr, vcm->group_j[0]);
 403             where();
 404             extract_binr(rb, icx, DIM*vcm->nr, vcm->group_x[0]);
 405             where();
 406             extract_binr(rb, ici, DIM*DIM*vcm->nr, vcm->group_i[0][0]);
 407             where();
 408         }
 409     }
 410
 411     if (nsig > 0)
 412     {
 413         extract_binr(rb, isig, nsig, sig);
 414     }
 415     where();
 416 }
 417
 418 int do_per_step(gmx_int64_t step, gmx_int64_t nstep)
 419 {
 420     if (nstep != 0)
 421     {
 422         return ((step % nstep) == 0);
 423     }
 424     else
 425     {
 426         return 0;
 427     }
 428 }