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