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 #if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
  41 /* _isnan() */
  42 #include <float.h>
  43 #endif
  44
  45 #include "typedefs.h"
  46 #include "string2.h"
  47 #include "smalloc.h"
  48 #include "mdrun.h"
  49 #include "domdec.h"
  50 #include "mtop_util.h"
  51 #include "gmx_wallcycle.h"
  52 #include "vcm.h"
  53 #include "nrnb.h"
  54
  55 /* Is the signal in one simulation independent of other simulations? */
  56 gmx_bool gs_simlocal[eglsNR] = { TRUE, FALSE, FALSE, TRUE };
  57
  58 /* check which of the multisim simulations has the shortest number of
  59    steps and return that number of nsteps */
  60 gmx_large_int_t get_multisim_nsteps(const t_commrec *cr,
  61                                     gmx_large_int_t nsteps)
  62 {
  63     gmx_large_int_t steps_out;
  64
  65     if MASTER(cr)
  66     {
  67         gmx_large_int_t *buf;
  68         int s;
  69
  70         snew(buf,cr->ms->nsim);
  71
  72         buf[cr->ms->sim] = nsteps;
  73         gmx_sumli_sim(cr->ms->nsim, buf, cr->ms);
  74
  75         steps_out=-1;
  76         for(s=0; s<cr->ms->nsim; s++)
  77         {
  78             /* find the smallest positive number */
  79             if (buf[s]>= 0 && ((steps_out < 0) || (buf[s]<steps_out)) )
  80             {
  81                 steps_out=buf[s];
  82             }
  83         }
  84         sfree(buf);
  85
  86         /* if we're the limiting simulation, don't do anything */
  87         if (steps_out>=0 && steps_out<nsteps)
  88         {
  89             char strbuf[255];
  90             snprintf(strbuf, 255, "Will stop simulation %%d after %s steps (another simulation will end then).\n", gmx_large_int_pfmt);
  91             fprintf(stderr, strbuf, cr->ms->sim, steps_out);
  92         }
  93     }
  94     /* broadcast to non-masters */
  95     gmx_bcast(sizeof(gmx_large_int_t), &steps_out, cr);
  96     return steps_out;
  97 }
  98
  99 int multisim_min(const gmx_multisim_t *ms,int nmin,int n)
 100 {
 101     int  *buf;
 102     gmx_bool bPos,bEqual;
 103     int  s,d;
 104
 105     snew(buf,ms->nsim);
 106     buf[ms->sim] = n;
 107     gmx_sumi_sim(ms->nsim,buf,ms);
 108     bPos   = TRUE;
 109     bEqual = TRUE;
 110     for(s=0; s<ms->nsim; s++)
 111     {
 112         bPos   = bPos   && (buf[s] > 0);
 113         bEqual = bEqual && (buf[s] == buf[0]);
 114     }
 115     if (bPos)
 116     {
 117         if (bEqual)
 118         {
 119             nmin = min(nmin,buf[0]);
 120         }
 121         else
 122         {
 123             /* Find the least common multiple */
 124             for(d=2; d<nmin; d++)
 125             {
 126                 s = 0;
 127                 while (s < ms->nsim && d % buf[s] == 0)
 128                 {
 129                     s++;
 130                 }
 131                 if (s == ms->nsim)
 132                 {
 133                     /* We found the LCM and it is less than nmin */
 134                     nmin = d;
 135                     break;
 136                 }
 137             }
 138         }
 139     }
 140     sfree(buf);
 141
 142     return nmin;
 143 }
 144
 145 int multisim_nstsimsync(const t_commrec *cr,
 146                         const t_inputrec *ir,int repl_ex_nst)
 147 {
 148     int nmin;
 149
 150     if (MASTER(cr))
 151     {
 152         nmin = INT_MAX;
 153         nmin = multisim_min(cr->ms,nmin,ir->nstlist);
 154         nmin = multisim_min(cr->ms,nmin,ir->nstcalcenergy);
 155         nmin = multisim_min(cr->ms,nmin,repl_ex_nst);
 156         if (nmin == INT_MAX)
 157         {
 158             gmx_fatal(FARGS,"Can not find an appropriate interval for inter-simulation communication, since nstlist, nstcalcenergy and -replex are all <= 0");
 159         }
 160         /* Avoid inter-simulation communication at every (second) step */
 161         if (nmin <= 2)
 162         {
 163             nmin = 10;
 164         }
 165     }
 166
 167     gmx_bcast(sizeof(int),&nmin,cr);
 168
 169     return nmin;
 170 }
 171
 172 void init_global_signals(globsig_t *gs,const t_commrec *cr,
 173                          const t_inputrec *ir,int repl_ex_nst)
 174 {
 175     int i;
 176
 177     if (MULTISIM(cr))
 178     {
 179         gs->nstms = multisim_nstsimsync(cr,ir,repl_ex_nst);
 180         if (debug)
 181         {
 182             fprintf(debug,"Syncing simulations for checkpointing and termination every %d steps\n",gs->nstms);
 183         }
 184     }
 185     else
 186     {
 187         gs->nstms = 1;
 188     }
 189
 190     for(i=0; i<eglsNR; i++)
 191     {
 192         gs->sig[i] = 0;
 193         gs->set[i] = 0;
 194     }
 195 }
 196
 197 void copy_coupling_state(t_state *statea,t_state *stateb,
 198                          gmx_ekindata_t *ekinda,gmx_ekindata_t *ekindb, t_grpopts* opts)
 199 {
 200
 201     /* MRS note -- might be able to get rid of some of the arguments.  Look over it when it's all debugged */
 202
 203     int i,j,nc;
 204
 205     /* Make sure we have enough space for x and v */
 206     if (statea->nalloc > stateb->nalloc)
 207     {
 208         stateb->nalloc = statea->nalloc;
 209         srenew(stateb->x,stateb->nalloc);
 210         srenew(stateb->v,stateb->nalloc);
 211     }
 212
 213     stateb->natoms     = statea->natoms;
 214     stateb->ngtc       = statea->ngtc;
 215     stateb->nnhpres    = statea->nnhpres;
 216     stateb->veta       = statea->veta;
 217     if (ekinda)
 218     {
 219         copy_mat(ekinda->ekin,ekindb->ekin);
 220         for (i=0; i<stateb->ngtc; i++)
 221         {
 222             ekindb->tcstat[i].T = ekinda->tcstat[i].T;
 223             ekindb->tcstat[i].Th = ekinda->tcstat[i].Th;
 224             copy_mat(ekinda->tcstat[i].ekinh,ekindb->tcstat[i].ekinh);
 225             copy_mat(ekinda->tcstat[i].ekinf,ekindb->tcstat[i].ekinf);
 226             ekindb->tcstat[i].ekinscalef_nhc =  ekinda->tcstat[i].ekinscalef_nhc;
 227             ekindb->tcstat[i].ekinscaleh_nhc =  ekinda->tcstat[i].ekinscaleh_nhc;
 228             ekindb->tcstat[i].vscale_nhc =  ekinda->tcstat[i].vscale_nhc;
 229         }
 230     }
 231     copy_rvecn(statea->x,stateb->x,0,stateb->natoms);
 232     copy_rvecn(statea->v,stateb->v,0,stateb->natoms);
 233     copy_mat(statea->box,stateb->box);
 234     copy_mat(statea->box_rel,stateb->box_rel);
 235     copy_mat(statea->boxv,stateb->boxv);
 236
 237     for (i = 0; i<stateb->ngtc; i++)
 238     {
 239         nc = i*opts->nhchainlength;
 240         for (j=0; j<opts->nhchainlength; j++)
 241         {
 242             stateb->nosehoover_xi[nc+j]  = statea->nosehoover_xi[nc+j];
 243             stateb->nosehoover_vxi[nc+j] = statea->nosehoover_vxi[nc+j];
 244         }
 245     }
 246     if (stateb->nhpres_xi != NULL)
 247     {
 248         for (i = 0; i<stateb->nnhpres; i++)
 249         {
 250             nc = i*opts->nhchainlength;
 251             for (j=0; j<opts->nhchainlength; j++)
 252             {
 253                 stateb->nhpres_xi[nc+j]  = statea->nhpres_xi[nc+j];
 254                 stateb->nhpres_vxi[nc+j] = statea->nhpres_vxi[nc+j];
 255             }
 256         }
 257     }
 258 }
 259
 260 real compute_conserved_from_auxiliary(t_inputrec *ir, t_state *state, t_extmass *MassQ)
 261 {
 262     real quantity = 0;
 263     switch (ir->etc)
 264     {
 265     case etcNO:
 266         break;
 267     case etcBERENDSEN:
 268         break;
 269     case etcNOSEHOOVER:
 270         quantity = NPT_energy(ir,state,MassQ);
 271         break;
 272     case etcVRESCALE:
 273         quantity = vrescale_energy(&(ir->opts),state->therm_integral);
 274         break;
 275     default:
 276         break;
 277     }
 278     return quantity;
 279 }
 280
 281 void compute_globals(FILE *fplog, gmx_global_stat_t gstat, t_commrec *cr, t_inputrec *ir,
 282                      t_forcerec *fr, gmx_ekindata_t *ekind,
 283                      t_state *state, t_state *state_global, t_mdatoms *mdatoms,
 284                      t_nrnb *nrnb, t_vcm *vcm, gmx_wallcycle_t wcycle,
 285                      gmx_enerdata_t *enerd,tensor force_vir, tensor shake_vir, tensor total_vir,
 286                      tensor pres, rvec mu_tot, gmx_constr_t constr,
 287                      globsig_t *gs,gmx_bool bInterSimGS,
 288                      matrix box, gmx_mtop_t *top_global, real *pcurr,
 289                      int natoms, gmx_bool *bSumEkinhOld, int flags)
 290 {
 291     int  i,gsi;
 292     real gs_buf[eglsNR];
 293     tensor corr_vir,corr_pres,shakeall_vir;
 294     gmx_bool bEner,bPres,bTemp, bVV;
 295     gmx_bool bRerunMD, bStopCM, bGStat, bIterate,
 296         bFirstIterate,bReadEkin,bEkinAveVel,bScaleEkin, bConstrain;
 297     real ekin,temp,prescorr,enercorr,dvdlcorr;
 298
 299     /* translate CGLO flags to gmx_booleans */
 300     bRerunMD = flags & CGLO_RERUNMD;
 301     bStopCM = flags & CGLO_STOPCM;
 302     bGStat = flags & CGLO_GSTAT;
 303
 304     bReadEkin = (flags & CGLO_READEKIN);
 305     bScaleEkin = (flags & CGLO_SCALEEKIN);
 306     bEner = flags & CGLO_ENERGY;
 307     bTemp = flags & CGLO_TEMPERATURE;
 308     bPres  = (flags & CGLO_PRESSURE);
 309     bConstrain = (flags & CGLO_CONSTRAINT);
 310     bIterate = (flags & CGLO_ITERATE);
 311     bFirstIterate = (flags & CGLO_FIRSTITERATE);
 312
 313     /* we calculate a full state kinetic energy either with full-step velocity verlet
 314        or half step where we need the pressure */
 315
 316     bEkinAveVel = (ir->eI==eiVV || (ir->eI==eiVVAK && bPres) || bReadEkin);
 317
 318     /* in initalization, it sums the shake virial in vv, and to
 319        sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
 320
 321     /* ########## Kinetic energy  ############## */
 322
 323     if (bTemp)
 324     {
 325         /* Non-equilibrium MD: this is parallellized, but only does communication
 326          * when there really is NEMD.
 327          */
 328
 329         if (PAR(cr) && (ekind->bNEMD))
 330         {
 331             accumulate_u(cr,&(ir->opts),ekind);
 332         }
 333         debug_gmx();
 334         if (bReadEkin)
 335         {
 336             restore_ekinstate_from_state(cr,ekind,&state_global->ekinstate);
 337         }
 338         else
 339         {
 340
 341             calc_ke_part(state,&(ir->opts),mdatoms,ekind,nrnb,bEkinAveVel,bIterate);
 342         }
 343
 344         debug_gmx();
 345
 346         /* Calculate center of mass velocity if necessary, also parallellized */
 347         if (bStopCM && !bRerunMD && bEner)
 348         {
 349             calc_vcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms,
 350                          state->x,state->v,vcm);
 351         }
 352     }
 353
 354     if (bTemp || bPres || bEner || bConstrain)
 355     {
 356         if (!bGStat)
 357         {
 358             /* We will not sum ekinh_old,
 359              * so signal that we still have to do it.
 360              */
 361             *bSumEkinhOld = TRUE;
 362
 363         }
 364         else
 365         {
 366             if (gs != NULL)
 367             {
 368                 for(i=0; i<eglsNR; i++)
 369                 {
 370                     gs_buf[i] = gs->sig[i];
 371                 }
 372             }
 373             if (PAR(cr))
 374             {
 375                 wallcycle_start(wcycle,ewcMoveE);
 376                 global_stat(fplog,gstat,cr,enerd,force_vir,shake_vir,mu_tot,
 377                             ir,ekind,constr,vcm,
 378                             gs != NULL ? eglsNR : 0,gs_buf,
 379                             top_global,state,
 380                             *bSumEkinhOld,flags);
 381                 wallcycle_stop(wcycle,ewcMoveE);
 382             }
 383             if (gs != NULL)
 384             {
 385                 if (MULTISIM(cr) && bInterSimGS)
 386                 {
 387                     if (MASTER(cr))
 388                     {
 389                         /* Communicate the signals between the simulations */
 390                         gmx_sum_sim(eglsNR,gs_buf,cr->ms);
 391                     }
 392                     /* Communicate the signals form the master to the others */
 393                     gmx_bcast(eglsNR*sizeof(gs_buf[0]),gs_buf,cr);
 394                 }
 395                 for(i=0; i<eglsNR; i++)
 396                 {
 397                     if (bInterSimGS || gs_simlocal[i])
 398                     {
 399                         /* Set the communicated signal only when it is non-zero,
 400                          * since signals might not be processed at each MD step.
 401                          */
 402                         gsi = (gs_buf[i] >= 0 ?
 403                                (int)(gs_buf[i] + 0.5) :
 404                                (int)(gs_buf[i] - 0.5));
 405                         if (gsi != 0)
 406                         {
 407                             gs->set[i] = gsi;
 408                         }
 409                         /* Turn off the local signal */
 410                         gs->sig[i] = 0;
 411                     }
 412                 }
 413             }
 414             *bSumEkinhOld = FALSE;
 415         }
 416     }
 417
 418     if (!ekind->bNEMD && debug && bTemp && (vcm->nr > 0))
 419     {
 420         correct_ekin(debug,
 421                      mdatoms->start,mdatoms->start+mdatoms->homenr,
 422                      state->v,vcm->group_p[0],
 423                      mdatoms->massT,mdatoms->tmass,ekind->ekin);
 424     }
 425
 426     if (bEner) {
 427         /* Do center of mass motion removal */
 428         if (bStopCM && !bRerunMD) /* is this correct?  Does it get called too often with this logic? */
 429         {
 430             check_cm_grp(fplog,vcm,ir,1);
 431             do_stopcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms->cVCM,
 432                           state->x,state->v,vcm);
 433             inc_nrnb(nrnb,eNR_STOPCM,mdatoms->homenr);
 434         }
 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,&(enerd->term[F_DKDL]),
 452                                        bEkinAveVel,bIterate,bScaleEkin);
 453
 454         enerd->term[F_EKIN] = trace(ekind->ekin);
 455     }
 456
 457     /* ##########  Long range energy information ###### */
 458
 459     if (bEner || bPres || bConstrain)
 460     {
 461         calc_dispcorr(fplog,ir,fr,0,top_global->natoms,box,state->lambda,
 462                       corr_pres,corr_vir,&prescorr,&enercorr,&dvdlcorr);
 463     }
 464
 465     if (bEner && bFirstIterate)
 466     {
 467         enerd->term[F_DISPCORR] = enercorr;
 468         enerd->term[F_EPOT] += enercorr;
 469         enerd->term[F_DVDL] += dvdlcorr;
 470         if (fr->efep != efepNO) {
 471             enerd->dvdl_lin += dvdlcorr;
 472         }
 473     }
 474
 475     /* ########## Now pressure ############## */
 476     if (bPres || bConstrain)
 477     {
 478
 479         m_add(force_vir,shake_vir,total_vir);
 480
 481         /* Calculate pressure and apply LR correction if PPPM is used.
 482          * Use the box from last timestep since we already called update().
 483          */
 484
 485         enerd->term[F_PRES] = calc_pres(fr->ePBC,ir->nwall,box,ekind->ekin,total_vir,pres,
 486                                         (fr->eeltype==eelPPPM)?enerd->term[F_COUL_RECIP]:0.0);
 487
 488         /* Calculate long range corrections to pressure and energy */
 489         /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
 490            and computes enerd->term[F_DISPCORR].  Also modifies the
 491            total_vir and pres tesors */
 492
 493         m_add(total_vir,corr_vir,total_vir);
 494         m_add(pres,corr_pres,pres);
 495         enerd->term[F_PDISPCORR] = prescorr;
 496         enerd->term[F_PRES] += prescorr;
 497         *pcurr = enerd->term[F_PRES];
 498         /* calculate temperature using virial */
 499         enerd->term[F_VTEMP] = calc_temp(trace(total_vir),ir->opts.nrdf[0]);
 500
 501     }
 502 }
 503
 504 void check_nst_param(FILE *fplog,t_commrec *cr,
 505                      const char *desc_nst,int nst,
 506                      const char *desc_p,int *p)
 507 {
 508     char buf[STRLEN];
 509
 510     if (*p > 0 && *p % nst != 0)
 511     {
 512         /* Round up to the next multiple of nst */
 513         *p = ((*p)/nst + 1)*nst;
 514         sprintf(buf,"NOTE: %s changes %s to %d\n",desc_nst,desc_p,*p);
 515         md_print_warning(cr,fplog,buf);
 516     }
 517 }
 518
 519 void reset_all_counters(FILE *fplog,t_commrec *cr,
 520                         gmx_large_int_t step,
 521                         gmx_large_int_t *step_rel,t_inputrec *ir,
 522                         gmx_wallcycle_t wcycle,t_nrnb *nrnb,
 523                         gmx_runtime_t *runtime)
 524 {
 525     char buf[STRLEN],sbuf[STEPSTRSIZE];
 526
 527     /* Reset all the counters related to performance over the run */
 528     sprintf(buf,"Step %s: resetting all time and cycle counters\n",
 529             gmx_step_str(step,sbuf));
 530     md_print_warning(cr,fplog,buf);
 531
 532     wallcycle_stop(wcycle,ewcRUN);
 533     wallcycle_reset_all(wcycle);
 534     if (DOMAINDECOMP(cr))
 535     {
 536         reset_dd_statistics_counters(cr->dd);
 537     }
 538     init_nrnb(nrnb);
 539     ir->init_step += *step_rel;
 540     ir->nsteps    -= *step_rel;
 541     *step_rel = 0;
 542     wallcycle_start(wcycle,ewcRUN);
 543     runtime_start(runtime);
 544     print_date_and_time(fplog,cr->nodeid,"Restarted time",runtime);
 545 }
 546
 547 void min_zero(int *n,int i)
 548 {
 549     if (i > 0 && (*n == 0 || i < *n))
 550     {
 551         *n = i;
 552     }
 553 }
 554
 555 int lcd4(int i1,int i2,int i3,int i4)
 556 {
 557     int nst;
 558
 559     nst = 0;
 560     min_zero(&nst,i1);
 561     min_zero(&nst,i2);
 562     min_zero(&nst,i3);
 563     min_zero(&nst,i4);
 564     if (nst == 0)
 565     {
 566         gmx_incons("All 4 inputs for determininig nstglobalcomm are <= 0");
 567     }
 568
 569     while (nst > 1 && ((i1 > 0 && i1 % nst != 0)  ||
 570                        (i2 > 0 && i2 % nst != 0)  ||
 571                        (i3 > 0 && i3 % nst != 0)  ||
 572                        (i4 > 0 && i4 % nst != 0)))
 573     {
 574         nst--;
 575     }
 576
 577     return nst;
 578 }
 579
 580 int check_nstglobalcomm(FILE *fplog,t_commrec *cr,
 581                         int nstglobalcomm,t_inputrec *ir)
 582 {
 583     char buf[STRLEN];
 584
 585     if (!EI_DYNAMICS(ir->eI))
 586     {
 587         nstglobalcomm = 1;
 588     }
 589
 590     if (nstglobalcomm == -1)
 591     {
 592         if (!(ir->nstcalcenergy > 0 ||
 593               ir->nstlist > 0 ||
 594               ir->etc != etcNO ||
 595               ir->epc != epcNO))
 596         {
 597             nstglobalcomm = 10;
 598             if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
 599             {
 600                 nstglobalcomm = ir->nstenergy;
 601             }
 602         }
 603         else
 604         {
 605             /* Ensure that we do timely global communication for
 606              * (possibly) each of the four following options.
 607              */
 608             nstglobalcomm = lcd4(ir->nstcalcenergy,
 609                                  ir->nstlist,
 610                                  ir->etc != etcNO ? ir->nsttcouple : 0,
 611                                  ir->epc != epcNO ? ir->nstpcouple : 0);
 612         }
 613     }
 614     else
 615     {
 616         if (ir->nstlist > 0 &&
 617             nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
 618         {
 619             nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
 620             sprintf(buf,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm);
 621             md_print_warning(cr,fplog,buf);
 622         }
 623         if (ir->nstcalcenergy > 0)
 624         {
 625             check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
 626                             "nstcalcenergy",&ir->nstcalcenergy);
 627         }
 628         if (ir->etc != etcNO && ir->nsttcouple > 0)
 629         {
 630             check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
 631                             "nsttcouple",&ir->nsttcouple);
 632         }
 633         if (ir->epc != epcNO && ir->nstpcouple > 0)
 634         {
 635             check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
 636                             "nstpcouple",&ir->nstpcouple);
 637         }
 638
 639         check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
 640                         "nstenergy",&ir->nstenergy);
 641
 642         check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
 643                         "nstlog",&ir->nstlog);
 644     }
 645
 646     if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
 647     {
 648         sprintf(buf,"WARNING: Changing nstcomm from %d to %d\n",
 649                 ir->nstcomm,nstglobalcomm);
 650         md_print_warning(cr,fplog,buf);
 651         ir->nstcomm = nstglobalcomm;
 652     }
 653
 654     return nstglobalcomm;
 655 }
 656
 657 void check_ir_old_tpx_versions(t_commrec *cr,FILE *fplog,
 658                                t_inputrec *ir,gmx_mtop_t *mtop)
 659 {
 660     /* Check required for old tpx files */
 661     if (IR_TWINRANGE(*ir) && ir->nstlist > 1 &&
 662         ir->nstcalcenergy % ir->nstlist != 0)
 663     {
 664         md_print_warning(cr,fplog,"Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies");
 665
 666         if (gmx_mtop_ftype_count(mtop,F_CONSTR) +
 667             gmx_mtop_ftype_count(mtop,F_CONSTRNC) > 0 &&
 668             ir->eConstrAlg == econtSHAKE)
 669         {
 670             md_print_warning(cr,fplog,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect");
 671             if (ir->epc != epcNO)
 672             {
 673                 gmx_fatal(FARGS,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
 674             }
 675         }
 676         check_nst_param(fplog,cr,"nstlist",ir->nstlist,
 677                         "nstcalcenergy",&ir->nstcalcenergy);
 678         if (ir->epc != epcNO)
 679         {
 680             check_nst_param(fplog,cr,"nstlist",ir->nstlist,
 681                             "nstpcouple",&ir->nstpcouple);
 682         }
 683         check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
 684                         "nstenergy",&ir->nstenergy);
 685         check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
 686                         "nstlog",&ir->nstlog);
 687         if (ir->efep != efepNO)
 688         {
 689             check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
 690                             "nstdhdl",&ir->nstdhdl);
 691         }
 692     }
 693 }
 694
 695 void rerun_parallel_comm(t_commrec *cr,t_trxframe *fr,
 696                          gmx_bool *bNotLastFrame)
 697 {
 698     gmx_bool bAlloc;
 699     rvec *xp,*vp;
 700
 701     bAlloc = (fr->natoms == 0);
 702
 703     if (MASTER(cr) && !*bNotLastFrame)
 704     {
 705         fr->natoms = -1;
 706     }
 707     xp = fr->x;
 708     vp = fr->v;
 709     gmx_bcast(sizeof(*fr),fr,cr);
 710     fr->x = xp;
 711     fr->v = vp;
 712
 713     *bNotLastFrame = (fr->natoms >= 0);
 714
 715     if (*bNotLastFrame && PARTDECOMP(cr))
 716     {
 717         /* x and v are the only variable size quantities stored in trr
 718          * that are required for rerun (f is not needed).
 719          */
 720         if (bAlloc)
 721         {
 722             snew(fr->x,fr->natoms);
 723             snew(fr->v,fr->natoms);
 724         }
 725         if (fr->bX)
 726         {
 727             gmx_bcast(fr->natoms*sizeof(fr->x[0]),fr->x[0],cr);
 728         }
 729         if (fr->bV)
 730         {
 731             gmx_bcast(fr->natoms*sizeof(fr->v[0]),fr->v[0],cr);
 732         }
 733     }
 734 }
 735
 736 void md_print_warning(const t_commrec *cr,FILE *fplog,const char *buf)
 737 {
 738     if (MASTER(cr))
 739     {
 740         fprintf(stderr,"\n%s\n",buf);
 741     }
 742     if (fplog)
 743     {
 744         fprintf(fplog,"\n%s\n",buf);
 745     }
 746 }