src/gromacs/mdlib/md_support.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include "typedefs.h"
  41 #include "string2.h"
  42 #include "smalloc.h"
  43 #include "mdrun.h"
  44 #include "domdec.h"
  45 #include "mtop_util.h"
  46 #include "vcm.h"
  47 #include "nrnb.h"
  48 #include "macros.h"
  49 #include "md_logging.h"
  50 #include "md_support.h"
  51
  52 #include "gromacs/timing/wallcycle.h"
  53
  54 /* Is the signal in one simulation independent of other simulations? */
  55 gmx_bool gs_simlocal[eglsNR] = { TRUE, FALSE, FALSE, TRUE };
  56
  57 /* check which of the multisim simulations has the shortest number of
  58    steps and return that number of nsteps */
  59 gmx_int64_t get_multisim_nsteps(const t_commrec *cr,
  60                                 gmx_int64_t      nsteps)
  61 {
  62     gmx_int64_t steps_out;
  63
  64     if (MASTER(cr))
  65     {
  66         gmx_int64_t     *buf;
  67         int              s;
  68
  69         snew(buf, cr->ms->nsim);
  70
  71         buf[cr->ms->sim] = nsteps;
  72         gmx_sumli_sim(cr->ms->nsim, buf, cr->ms);
  73
  74         steps_out = -1;
  75         for (s = 0; s < cr->ms->nsim; s++)
  76         {
  77             /* find the smallest positive number */
  78             if (buf[s] >= 0 && ((steps_out < 0) || (buf[s] < steps_out)) )
  79             {
  80                 steps_out = buf[s];
  81             }
  82         }
  83         sfree(buf);
  84
  85         /* if we're the limiting simulation, don't do anything */
  86         if (steps_out >= 0 && steps_out < nsteps)
  87         {
  88             char strbuf[255];
  89             snprintf(strbuf, 255, "Will stop simulation %%d after %s steps (another simulation will end then).\n", "%"GMX_PRId64);
  90             fprintf(stderr, strbuf, cr->ms->sim, steps_out);
  91         }
  92     }
  93     /* broadcast to non-masters */
  94     gmx_bcast(sizeof(gmx_int64_t), &steps_out, cr);
  95     return steps_out;
  96 }
  97
  98 int multisim_min(const gmx_multisim_t *ms, int nmin, int n)
  99 {
 100     int     *buf;
 101     gmx_bool bPos, bEqual;
 102     int      s, d;
 103
 104     snew(buf, ms->nsim);
 105     buf[ms->sim] = n;
 106     gmx_sumi_sim(ms->nsim, buf, ms);
 107     bPos   = TRUE;
 108     bEqual = TRUE;
 109     for (s = 0; s < ms->nsim; s++)
 110     {
 111         bPos   = bPos   && (buf[s] > 0);
 112         bEqual = bEqual && (buf[s] == buf[0]);
 113     }
 114     if (bPos)
 115     {
 116         if (bEqual)
 117         {
 118             nmin = min(nmin, buf[0]);
 119         }
 120         else
 121         {
 122             /* Find the least common multiple */
 123             for (d = 2; d < nmin; d++)
 124             {
 125                 s = 0;
 126                 while (s < ms->nsim && d % buf[s] == 0)
 127                 {
 128                     s++;
 129                 }
 130                 if (s == ms->nsim)
 131                 {
 132                     /* We found the LCM and it is less than nmin */
 133                     nmin = d;
 134                     break;
 135                 }
 136             }
 137         }
 138     }
 139     sfree(buf);
 140
 141     return nmin;
 142 }
 143
 144 int multisim_nstsimsync(const t_commrec *cr,
 145                         const t_inputrec *ir, int repl_ex_nst)
 146 {
 147     int nmin;
 148
 149     if (MASTER(cr))
 150     {
 151         nmin = INT_MAX;
 152         nmin = multisim_min(cr->ms, nmin, ir->nstlist);
 153         nmin = multisim_min(cr->ms, nmin, ir->nstcalcenergy);
 154         nmin = multisim_min(cr->ms, nmin, repl_ex_nst);
 155         if (nmin == INT_MAX)
 156         {
 157             gmx_fatal(FARGS, "Can not find an appropriate interval for inter-simulation communication, since nstlist, nstcalcenergy and -replex are all <= 0");
 158         }
 159         /* Avoid inter-simulation communication at every (second) step */
 160         if (nmin <= 2)
 161         {
 162             nmin = 10;
 163         }
 164     }
 165
 166     gmx_bcast(sizeof(int), &nmin, cr);
 167
 168     return nmin;
 169 }
 170
 171 void init_global_signals(globsig_t *gs, const t_commrec *cr,
 172                          const t_inputrec *ir, int repl_ex_nst)
 173 {
 174     int i;
 175
 176     if (MULTISIM(cr))
 177     {
 178         gs->nstms = multisim_nstsimsync(cr, ir, repl_ex_nst);
 179         if (debug)
 180         {
 181             fprintf(debug, "Syncing simulations for checkpointing and termination every %d steps\n", gs->nstms);
 182         }
 183     }
 184     else
 185     {
 186         gs->nstms = 1;
 187     }
 188
 189     for (i = 0; i < eglsNR; i++)
 190     {
 191         gs->sig[i] = 0;
 192         gs->set[i] = 0;
 193     }
 194 }
 195
 196 void copy_coupling_state(t_state *statea, t_state *stateb,
 197                          gmx_ekindata_t *ekinda, gmx_ekindata_t *ekindb, t_grpopts* opts)
 198 {
 199
 200     /* MRS note -- might be able to get rid of some of the arguments.  Look over it when it's all debugged */
 201
 202     int i, j, nc;
 203
 204     /* Make sure we have enough space for x and v */
 205     if (statea->nalloc > stateb->nalloc)
 206     {
 207         stateb->nalloc = statea->nalloc;
 208         srenew(stateb->x, stateb->nalloc);
 209         srenew(stateb->v, stateb->nalloc);
 210     }
 211
 212     stateb->natoms     = statea->natoms;
 213     stateb->ngtc       = statea->ngtc;
 214     stateb->nnhpres    = statea->nnhpres;
 215     stateb->veta       = statea->veta;
 216     if (ekinda)
 217     {
 218         copy_mat(ekinda->ekin, ekindb->ekin);
 219         for (i = 0; i < stateb->ngtc; i++)
 220         {
 221             ekindb->tcstat[i].T  = ekinda->tcstat[i].T;
 222             ekindb->tcstat[i].Th = ekinda->tcstat[i].Th;
 223             copy_mat(ekinda->tcstat[i].ekinh, ekindb->tcstat[i].ekinh);
 224             copy_mat(ekinda->tcstat[i].ekinf, ekindb->tcstat[i].ekinf);
 225             ekindb->tcstat[i].ekinscalef_nhc =  ekinda->tcstat[i].ekinscalef_nhc;
 226             ekindb->tcstat[i].ekinscaleh_nhc =  ekinda->tcstat[i].ekinscaleh_nhc;
 227             ekindb->tcstat[i].vscale_nhc     =  ekinda->tcstat[i].vscale_nhc;
 228         }
 229     }
 230     copy_rvecn(statea->x, stateb->x, 0, stateb->natoms);
 231     copy_rvecn(statea->v, stateb->v, 0, stateb->natoms);
 232     copy_mat(statea->box, stateb->box);
 233     copy_mat(statea->box_rel, stateb->box_rel);
 234     copy_mat(statea->boxv, stateb->boxv);
 235
 236     for (i = 0; i < stateb->ngtc; i++)
 237     {
 238         nc = i*opts->nhchainlength;
 239         for (j = 0; j < opts->nhchainlength; j++)
 240         {
 241             stateb->nosehoover_xi[nc+j]  = statea->nosehoover_xi[nc+j];
 242             stateb->nosehoover_vxi[nc+j] = statea->nosehoover_vxi[nc+j];
 243         }
 244     }
 245     if (stateb->nhpres_xi != NULL)
 246     {
 247         for (i = 0; i < stateb->nnhpres; i++)
 248         {
 249             nc = i*opts->nhchainlength;
 250             for (j = 0; j < opts->nhchainlength; j++)
 251             {
 252                 stateb->nhpres_xi[nc+j]  = statea->nhpres_xi[nc+j];
 253                 stateb->nhpres_vxi[nc+j] = statea->nhpres_vxi[nc+j];
 254             }
 255         }
 256     }
 257 }
 258
 259 real compute_conserved_from_auxiliary(t_inputrec *ir, t_state *state, t_extmass *MassQ)
 260 {
 261     real quantity = 0;
 262     switch (ir->etc)
 263     {
 264         case etcNO:
 265             break;
 266         case etcBERENDSEN:
 267             break;
 268         case etcNOSEHOOVER:
 269             quantity = NPT_energy(ir, state, MassQ);
 270             break;
 271         case etcVRESCALE:
 272             quantity = vrescale_energy(&(ir->opts), state->therm_integral);
 273             break;
 274         default:
 275             break;
 276     }
 277     return quantity;
 278 }
 279
 280 void compute_globals(FILE *fplog, gmx_global_stat_t gstat, t_commrec *cr, t_inputrec *ir,
 281                      t_forcerec *fr, gmx_ekindata_t *ekind,
 282                      t_state *state, t_state *state_global, t_mdatoms *mdatoms,
 283                      t_nrnb *nrnb, t_vcm *vcm, gmx_wallcycle_t wcycle,
 284                      gmx_enerdata_t *enerd, tensor force_vir, tensor shake_vir, tensor total_vir,
 285                      tensor pres, rvec mu_tot, gmx_constr_t constr,
 286                      globsig_t *gs, gmx_bool bInterSimGS,
 287                      matrix box, gmx_mtop_t *top_global,
 288                      gmx_bool *bSumEkinhOld, int flags)
 289 {
 290     int      i, gsi;
 291     real     gs_buf[eglsNR];
 292     tensor   corr_vir, corr_pres;
 293     gmx_bool bEner, bPres, bTemp, bVV;
 294     gmx_bool bRerunMD, bStopCM, bGStat, bIterate,
 295              bFirstIterate, bReadEkin, bEkinAveVel, bScaleEkin, bConstrain;
 296     real     ekin, temp, prescorr, enercorr, dvdlcorr, dvdl_ekin;
 297
 298     /* translate CGLO flags to gmx_booleans */
 299     bRerunMD = flags & CGLO_RERUNMD;
 300     bStopCM  = flags & CGLO_STOPCM;
 301     bGStat   = flags & CGLO_GSTAT;
 302
 303     bReadEkin     = (flags & CGLO_READEKIN);
 304     bScaleEkin    = (flags & CGLO_SCALEEKIN);
 305     bEner         = flags & CGLO_ENERGY;
 306     bTemp         = flags & CGLO_TEMPERATURE;
 307     bPres         = (flags & CGLO_PRESSURE);
 308     bConstrain    = (flags & CGLO_CONSTRAINT);
 309     bIterate      = (flags & CGLO_ITERATE);
 310     bFirstIterate = (flags & CGLO_FIRSTITERATE);
 311
 312     /* we calculate a full state kinetic energy either with full-step velocity verlet
 313        or half step where we need the pressure */
 314
 315     bEkinAveVel = (ir->eI == eiVV || (ir->eI == eiVVAK && bPres) || bReadEkin);
 316
 317     /* in initalization, it sums the shake virial in vv, and to
 318        sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
 319
 320     /* ########## Kinetic energy  ############## */
 321
 322     if (bTemp)
 323     {
 324         /* Non-equilibrium MD: this is parallellized, but only does communication
 325          * when there really is NEMD.
 326          */
 327
 328         if (PAR(cr) && (ekind->bNEMD))
 329         {
 330             accumulate_u(cr, &(ir->opts), ekind);
 331         }
 332         debug_gmx();
 333         if (bReadEkin)
 334         {
 335             restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
 336         }
 337         else
 338         {
 339
 340             calc_ke_part(state, &(ir->opts), mdatoms, ekind, nrnb, bEkinAveVel, bIterate);
 341         }
 342
 343         debug_gmx();
 344     }
 345
 346     /* Calculate center of mass velocity if necessary, also parallellized */
 347     if (bStopCM)
 348     {
 349         calc_vcm_grp(mdatoms->start, mdatoms->homenr, mdatoms,
 350                      state->x, state->v, vcm);
 351     }
 352
 353     if (bTemp || bStopCM || bPres || bEner || bConstrain)
 354     {
 355         if (!bGStat)
 356         {
 357             /* We will not sum ekinh_old,
 358              * so signal that we still have to do it.
 359              */
 360             *bSumEkinhOld = TRUE;
 361
 362         }
 363         else
 364         {
 365             if (gs != NULL)
 366             {
 367                 for (i = 0; i < eglsNR; i++)
 368                 {
 369                     gs_buf[i] = gs->sig[i];
 370                 }
 371             }
 372             if (PAR(cr))
 373             {
 374                 wallcycle_start(wcycle, ewcMoveE);
 375                 global_stat(fplog, gstat, cr, enerd, force_vir, shake_vir, mu_tot,
 376                             ir, ekind, constr, bStopCM ? vcm : NULL,
 377                             gs != NULL ? eglsNR : 0, gs_buf,
 378                             top_global, state,
 379                             *bSumEkinhOld, flags);
 380                 wallcycle_stop(wcycle, ewcMoveE);
 381             }
 382             if (gs != NULL)
 383             {
 384                 if (MULTISIM(cr) && bInterSimGS)
 385                 {
 386                     if (MASTER(cr))
 387                     {
 388                         /* Communicate the signals between the simulations */
 389                         gmx_sum_sim(eglsNR, gs_buf, cr->ms);
 390                     }
 391                     /* Communicate the signals form the master to the others */
 392                     gmx_bcast(eglsNR*sizeof(gs_buf[0]), gs_buf, cr);
 393                 }
 394                 for (i = 0; i < eglsNR; i++)
 395                 {
 396                     if (bInterSimGS || gs_simlocal[i])
 397                     {
 398                         /* Set the communicated signal only when it is non-zero,
 399                          * since signals might not be processed at each MD step.
 400                          */
 401                         gsi = (gs_buf[i] >= 0 ?
 402                                (int)(gs_buf[i] + 0.5) :
 403                                (int)(gs_buf[i] - 0.5));
 404                         if (gsi != 0)
 405                         {
 406                             gs->set[i] = gsi;
 407                         }
 408                         /* Turn off the local signal */
 409                         gs->sig[i] = 0;
 410                     }
 411                 }
 412             }
 413             *bSumEkinhOld = FALSE;
 414         }
 415     }
 416
 417     if (!ekind->bNEMD && debug && bTemp && (vcm->nr > 0))
 418     {
 419         correct_ekin(debug,
 420                      mdatoms->start, mdatoms->start+mdatoms->homenr,
 421                      state->v, vcm->group_p[0],
 422                      mdatoms->massT, mdatoms->tmass, ekind->ekin);
 423     }
 424
 425     /* Do center of mass motion removal */
 426     if (bStopCM)
 427     {
 428         check_cm_grp(fplog, vcm, ir, 1);
 429         do_stopcm_grp(mdatoms->start, mdatoms->homenr, mdatoms->cVCM,
 430                       state->x, state->v, vcm);
 431         inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
 432     }
 433
 434     if (bEner)
 435     {
 436         /* Calculate the amplitude of the cosine velocity profile */
 437         ekind->cosacc.vcos = ekind->cosacc.mvcos/mdatoms->tmass;
 438     }
 439
 440     if (bTemp)
 441     {
 442         /* Sum the kinetic energies of the groups & calc temp */
 443         /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with IR_NPT_TROTTER */
 444         /* three maincase:  VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
 445            Leap with AveVel is not supported; it's not clear that it will actually work.
 446            bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
 447            If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
 448            bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
 449            If FALSE, we go ahead and erase over it.
 450          */
 451         enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, &dvdl_ekin,
 452                                        bEkinAveVel, bScaleEkin);
 453         enerd->dvdl_lin[efptMASS] = (double) dvdl_ekin;
 454
 455         enerd->term[F_EKIN] = trace(ekind->ekin);
 456     }
 457
 458     /* ##########  Long range energy information ###### */
 459
 460     if (bEner || bPres || bConstrain)
 461     {
 462         calc_dispcorr(fplog, ir, fr, 0, top_global->natoms, box, state->lambda[efptVDW],
 463                       corr_pres, corr_vir, &prescorr, &enercorr, &dvdlcorr);
 464     }
 465
 466     if (bEner && bFirstIterate)
 467     {
 468         enerd->term[F_DISPCORR]  = enercorr;
 469         enerd->term[F_EPOT]     += enercorr;
 470         enerd->term[F_DVDL_VDW] += dvdlcorr;
 471     }
 472
 473     /* ########## Now pressure ############## */
 474     if (bPres || bConstrain)
 475     {
 476
 477         m_add(force_vir, shake_vir, total_vir);
 478
 479         /* Calculate pressure and apply LR correction if PPPM is used.
 480          * Use the box from last timestep since we already called update().
 481          */
 482
 483         enerd->term[F_PRES] = calc_pres(fr->ePBC, ir->nwall, box, ekind->ekin, total_vir, pres);
 484
 485         /* Calculate long range corrections to pressure and energy */
 486         /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
 487            and computes enerd->term[F_DISPCORR].  Also modifies the
 488            total_vir and pres tesors */
 489
 490         m_add(total_vir, corr_vir, total_vir);
 491         m_add(pres, corr_pres, pres);
 492         enerd->term[F_PDISPCORR] = prescorr;
 493         enerd->term[F_PRES]     += prescorr;
 494     }
 495 }
 496
 497 void check_nst_param(FILE *fplog, t_commrec *cr,
 498                      const char *desc_nst, int nst,
 499                      const char *desc_p, int *p)
 500 {
 501     if (*p > 0 && *p % nst != 0)
 502     {
 503         /* Round up to the next multiple of nst */
 504         *p = ((*p)/nst + 1)*nst;
 505         md_print_warn(cr, fplog,
 506                       "NOTE: %s changes %s to %d\n", desc_nst, desc_p, *p);
 507     }
 508 }
 509
 510 void set_current_lambdas(gmx_int64_t step, t_lambda *fepvals, gmx_bool bRerunMD,
 511                          t_trxframe *rerun_fr, t_state *state_global, t_state *state, double lam0[])
 512 /* find the current lambdas.  If rerunning, we either read in a state, or a lambda value,
 513    requiring different logic. */
 514 {
 515     real frac;
 516     int  i, fep_state = 0;
 517     if (bRerunMD)
 518     {
 519         if (rerun_fr->bLambda)
 520         {
 521             if (fepvals->delta_lambda == 0)
 522             {
 523                 state_global->lambda[efptFEP] = rerun_fr->lambda;
 524                 for (i = 0; i < efptNR; i++)
 525                 {
 526                     if (i != efptFEP)
 527                     {
 528                         state->lambda[i] = state_global->lambda[i];
 529                     }
 530                 }
 531             }
 532             else
 533             {
 534                 /* find out between which two value of lambda we should be */
 535                 frac      = (step*fepvals->delta_lambda);
 536                 fep_state = floor(frac*fepvals->n_lambda);
 537                 /* interpolate between this state and the next */
 538                 /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
 539                 frac = (frac*fepvals->n_lambda)-fep_state;
 540                 for (i = 0; i < efptNR; i++)
 541                 {
 542                     state_global->lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state]) +
 543                         frac*(fepvals->all_lambda[i][fep_state+1]-fepvals->all_lambda[i][fep_state]);
 544                 }
 545             }
 546         }
 547         else if (rerun_fr->bFepState)
 548         {
 549             state_global->fep_state = rerun_fr->fep_state;
 550             for (i = 0; i < efptNR; i++)
 551             {
 552                 state_global->lambda[i] = fepvals->all_lambda[i][fep_state];
 553             }
 554         }
 555     }
 556     else
 557     {
 558         if (fepvals->delta_lambda != 0)
 559         {
 560             /* find out between which two value of lambda we should be */
 561             frac = (step*fepvals->delta_lambda);
 562             if (fepvals->n_lambda > 0)
 563             {
 564                 fep_state = floor(frac*fepvals->n_lambda);
 565                 /* interpolate between this state and the next */
 566                 /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
 567                 frac = (frac*fepvals->n_lambda)-fep_state;
 568                 for (i = 0; i < efptNR; i++)
 569                 {
 570                     state_global->lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state]) +
 571                         frac*(fepvals->all_lambda[i][fep_state+1]-fepvals->all_lambda[i][fep_state]);
 572                 }
 573             }
 574             else
 575             {
 576                 for (i = 0; i < efptNR; i++)
 577                 {
 578                     state_global->lambda[i] = lam0[i] + frac;
 579                 }
 580             }
 581         }
 582         else
 583         {
 584             if (state->fep_state > 0)
 585             {
 586                 state_global->fep_state = state->fep_state; /* state->fep is the one updated by bExpanded */
 587                 for (i = 0; i < efptNR; i++)
 588                 {
 589                     state_global->lambda[i] = fepvals->all_lambda[i][state_global->fep_state];
 590                 }
 591             }
 592         }
 593     }
 594     for (i = 0; i < efptNR; i++)
 595     {
 596         state->lambda[i] = state_global->lambda[i];
 597     }
 598 }
 599
 600 static void min_zero(int *n, int i)
 601 {
 602     if (i > 0 && (*n == 0 || i < *n))
 603     {
 604         *n = i;
 605     }
 606 }
 607
 608 static int lcd4(int i1, int i2, int i3, int i4)
 609 {
 610     int nst;
 611
 612     nst = 0;
 613     min_zero(&nst, i1);
 614     min_zero(&nst, i2);
 615     min_zero(&nst, i3);
 616     min_zero(&nst, i4);
 617     if (nst == 0)
 618     {
 619         gmx_incons("All 4 inputs for determininig nstglobalcomm are <= 0");
 620     }
 621
 622     while (nst > 1 && ((i1 > 0 && i1 % nst != 0)  ||
 623                        (i2 > 0 && i2 % nst != 0)  ||
 624                        (i3 > 0 && i3 % nst != 0)  ||
 625                        (i4 > 0 && i4 % nst != 0)))
 626     {
 627         nst--;
 628     }
 629
 630     return nst;
 631 }
 632
 633 int check_nstglobalcomm(FILE *fplog, t_commrec *cr,
 634                         int nstglobalcomm, t_inputrec *ir)
 635 {
 636     if (!EI_DYNAMICS(ir->eI))
 637     {
 638         nstglobalcomm = 1;
 639     }
 640
 641     if (nstglobalcomm == -1)
 642     {
 643         if (!(ir->nstcalcenergy > 0 ||
 644               ir->nstlist > 0 ||
 645               ir->etc != etcNO ||
 646               ir->epc != epcNO))
 647         {
 648             nstglobalcomm = 10;
 649             if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
 650             {
 651                 nstglobalcomm = ir->nstenergy;
 652             }
 653         }
 654         else
 655         {
 656             /* Ensure that we do timely global communication for
 657              * (possibly) each of the four following options.
 658              */
 659             nstglobalcomm = lcd4(ir->nstcalcenergy,
 660                                  ir->nstlist,
 661                                  ir->etc != etcNO ? ir->nsttcouple : 0,
 662                                  ir->epc != epcNO ? ir->nstpcouple : 0);
 663         }
 664     }
 665     else
 666     {
 667         if (ir->nstlist > 0 &&
 668             nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
 669         {
 670             nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
 671             md_print_warn(cr, fplog, "WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n", nstglobalcomm);
 672         }
 673         if (ir->nstcalcenergy > 0)
 674         {
 675             check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 676                             "nstcalcenergy", &ir->nstcalcenergy);
 677         }
 678         if (ir->etc != etcNO && ir->nsttcouple > 0)
 679         {
 680             check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 681                             "nsttcouple", &ir->nsttcouple);
 682         }
 683         if (ir->epc != epcNO && ir->nstpcouple > 0)
 684         {
 685             check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 686                             "nstpcouple", &ir->nstpcouple);
 687         }
 688
 689         check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 690                         "nstenergy", &ir->nstenergy);
 691
 692         check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 693                         "nstlog", &ir->nstlog);
 694     }
 695
 696     if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
 697     {
 698         md_print_warn(cr, fplog, "WARNING: Changing nstcomm from %d to %d\n",
 699                       ir->nstcomm, nstglobalcomm);
 700         ir->nstcomm = nstglobalcomm;
 701     }
 702
 703     return nstglobalcomm;
 704 }
 705
 706 void check_ir_old_tpx_versions(t_commrec *cr, FILE *fplog,
 707                                t_inputrec *ir, gmx_mtop_t *mtop)
 708 {
 709     /* Check required for old tpx files */
 710     if (IR_TWINRANGE(*ir) && ir->nstlist > 1 &&
 711         ir->nstcalcenergy % ir->nstlist != 0)
 712     {
 713         md_print_warn(cr, fplog, "Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies\n");
 714
 715         if (gmx_mtop_ftype_count(mtop, F_CONSTR) +
 716             gmx_mtop_ftype_count(mtop, F_CONSTRNC) > 0 &&
 717             ir->eConstrAlg == econtSHAKE)
 718         {
 719             md_print_warn(cr, fplog, "With twin-range cut-off's and SHAKE the virial and pressure are incorrect\n");
 720             if (ir->epc != epcNO)
 721             {
 722                 gmx_fatal(FARGS, "Can not do pressure coupling with twin-range cut-off's and SHAKE");
 723             }
 724         }
 725         check_nst_param(fplog, cr, "nstlist", ir->nstlist,
 726                         "nstcalcenergy", &ir->nstcalcenergy);
 727         if (ir->epc != epcNO)
 728         {
 729             check_nst_param(fplog, cr, "nstlist", ir->nstlist,
 730                             "nstpcouple", &ir->nstpcouple);
 731         }
 732         check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
 733                         "nstenergy", &ir->nstenergy);
 734         check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
 735                         "nstlog", &ir->nstlog);
 736         if (ir->efep != efepNO)
 737         {
 738             check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
 739                             "nstdhdl", &ir->fepvals->nstdhdl);
 740         }
 741     }
 742 }
 743
 744 void rerun_parallel_comm(t_commrec *cr, t_trxframe *fr,
 745                          gmx_bool *bNotLastFrame)
 746 {
 747     gmx_bool bAlloc;
 748     rvec    *xp, *vp;
 749
 750     bAlloc = (fr->natoms == 0);
 751
 752     if (MASTER(cr) && !*bNotLastFrame)
 753     {
 754         fr->natoms = -1;
 755     }
 756     xp = fr->x;
 757     vp = fr->v;
 758     gmx_bcast(sizeof(*fr), fr, cr);
 759     fr->x = xp;
 760     fr->v = vp;
 761
 762     *bNotLastFrame = (fr->natoms >= 0);
 763
 764     if (*bNotLastFrame && PARTDECOMP(cr))
 765     {
 766         /* x and v are the only variable size quantities stored in trr
 767          * that are required for rerun (f is not needed).
 768          */
 769         if (bAlloc)
 770         {
 771             snew(fr->x, fr->natoms);
 772             snew(fr->v, fr->natoms);
 773         }
 774         if (fr->bX)
 775         {
 776             gmx_bcast(fr->natoms*sizeof(fr->x[0]), fr->x[0], cr);
 777         }
 778         if (fr->bV)
 779         {
 780             gmx_bcast(fr->natoms*sizeof(fr->v[0]), fr->v[0], cr);
 781         }
 782     }
 783 }