Bug Summary

File:gromacs/mdlib/md_support.c
Location:line 301, column 5
Description:Value stored to 'bRerunMD' is never read

Annotated Source Code

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