1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
31 * For more info, check our website at http://www.gromacs.org
34 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
43 #if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
80 #include "mpelogging.h"
87 #include "compute_io.h"
89 #include "checkpoint.h"
90 #include "mtop_util.h"
91 #include "sighandler.h"
101 #include "corewrap.h"
106 enum { eglsNABNSB
, eglsCHKPT
, eglsTERM
, eglsRESETCOUNTERS
, eglsNR
};
107 /* Is the signal in one simulation independent of other simulations? */
108 bool gs_simlocal
[eglsNR
] = { TRUE
, FALSE
, FALSE
, TRUE
};
111 int nstms
; /* The frequency for intersimulation communication */
112 int sig
[eglsNR
]; /* The signal set by one process in do_md */
113 int set
[eglsNR
]; /* The communicated signal, equal for all processes */
117 static int multisim_min(const gmx_multisim_t
*ms
,int nmin
,int n
)
125 gmx_sumi_sim(ms
->nsim
,buf
,ms
);
128 for(s
=0; s
<ms
->nsim
; s
++)
130 bPos
= bPos
&& (buf
[s
] > 0);
131 bEqual
= bEqual
&& (buf
[s
] == buf
[0]);
137 nmin
= min(nmin
,buf
[0]);
141 /* Find the least common multiple */
142 for(d
=2; d
<nmin
; d
++)
145 while (s
< ms
->nsim
&& d
% buf
[s
] == 0)
151 /* We found the LCM and it is less than nmin */
163 static int multisim_nstsimsync(const t_commrec
*cr
,
164 const t_inputrec
*ir
,int repl_ex_nst
)
171 nmin
= multisim_min(cr
->ms
,nmin
,ir
->nstlist
);
172 nmin
= multisim_min(cr
->ms
,nmin
,ir
->nstcalcenergy
);
173 nmin
= multisim_min(cr
->ms
,nmin
,repl_ex_nst
);
176 gmx_fatal(FARGS
,"Can not find an appropriate interval for inter-simulation communication, since nstlist, nstcalcenergy and -replex are all <= 0");
178 /* Avoid inter-simulation communication at every (second) step */
185 gmx_bcast(sizeof(int),&nmin
,cr
);
190 static void init_global_signals(globsig_t
*gs
,const t_commrec
*cr
,
191 const t_inputrec
*ir
,int repl_ex_nst
)
197 gs
->nstms
= multisim_nstsimsync(cr
,ir
,repl_ex_nst
);
200 fprintf(debug
,"Syncing simulations for checkpointing and termination every %d steps\n",gs
->nstms
);
208 for(i
=0; i
<eglsNR
; i
++)
215 static void copy_coupling_state(t_state
*statea
,t_state
*stateb
,
216 gmx_ekindata_t
*ekinda
,gmx_ekindata_t
*ekindb
, t_grpopts
* opts
)
219 /* MRS note -- might be able to get rid of some of the arguments. Look over it when it's all debugged */
223 /* Make sure we have enough space for x and v */
224 if (statea
->nalloc
> stateb
->nalloc
)
226 stateb
->nalloc
= statea
->nalloc
;
227 srenew(stateb
->x
,stateb
->nalloc
);
228 srenew(stateb
->v
,stateb
->nalloc
);
231 stateb
->natoms
= statea
->natoms
;
232 stateb
->ngtc
= statea
->ngtc
;
233 stateb
->nnhpres
= statea
->nnhpres
;
234 stateb
->veta
= statea
->veta
;
237 copy_mat(ekinda
->ekin
,ekindb
->ekin
);
238 for (i
=0; i
<stateb
->ngtc
; i
++)
240 ekindb
->tcstat
[i
].T
= ekinda
->tcstat
[i
].T
;
241 ekindb
->tcstat
[i
].Th
= ekinda
->tcstat
[i
].Th
;
242 copy_mat(ekinda
->tcstat
[i
].ekinh
,ekindb
->tcstat
[i
].ekinh
);
243 copy_mat(ekinda
->tcstat
[i
].ekinf
,ekindb
->tcstat
[i
].ekinf
);
244 ekindb
->tcstat
[i
].ekinscalef_nhc
= ekinda
->tcstat
[i
].ekinscalef_nhc
;
245 ekindb
->tcstat
[i
].ekinscaleh_nhc
= ekinda
->tcstat
[i
].ekinscaleh_nhc
;
246 ekindb
->tcstat
[i
].vscale_nhc
= ekinda
->tcstat
[i
].vscale_nhc
;
249 copy_rvecn(statea
->x
,stateb
->x
,0,stateb
->natoms
);
250 copy_rvecn(statea
->v
,stateb
->v
,0,stateb
->natoms
);
251 copy_mat(statea
->box
,stateb
->box
);
252 copy_mat(statea
->box_rel
,stateb
->box_rel
);
253 copy_mat(statea
->boxv
,stateb
->boxv
);
255 for (i
= 0; i
<stateb
->ngtc
; i
++)
257 nc
= i
*opts
->nhchainlength
;
258 for (j
=0; j
<opts
->nhchainlength
; j
++)
260 stateb
->nosehoover_xi
[nc
+j
] = statea
->nosehoover_xi
[nc
+j
];
261 stateb
->nosehoover_vxi
[nc
+j
] = statea
->nosehoover_vxi
[nc
+j
];
264 if (stateb
->nhpres_xi
!= NULL
)
266 for (i
= 0; i
<stateb
->nnhpres
; i
++)
268 nc
= i
*opts
->nhchainlength
;
269 for (j
=0; j
<opts
->nhchainlength
; j
++)
271 stateb
->nhpres_xi
[nc
+j
] = statea
->nhpres_xi
[nc
+j
];
272 stateb
->nhpres_vxi
[nc
+j
] = statea
->nhpres_vxi
[nc
+j
];
278 static void compute_globals(FILE *fplog
, gmx_global_stat_t gstat
, t_commrec
*cr
, t_inputrec
*ir
,
279 t_forcerec
*fr
, gmx_ekindata_t
*ekind
,
280 t_state
*state
, t_state
*state_global
, t_mdatoms
*mdatoms
,
281 t_nrnb
*nrnb
, t_vcm
*vcm
, gmx_wallcycle_t wcycle
,
282 gmx_enerdata_t
*enerd
,tensor force_vir
, tensor shake_vir
, tensor total_vir
,
283 tensor pres
, rvec mu_tot
, gmx_constr_t constr
,
284 globsig_t
*gs
,bool bInterSimGS
,
285 matrix box
, gmx_mtop_t
*top_global
, real
*pcurr
,
286 int natoms
, bool *bSumEkinhOld
, int flags
)
290 tensor corr_vir
,corr_pres
,shakeall_vir
;
291 bool bEner
,bPres
,bTemp
, bVV
;
292 bool bRerunMD
, bStopCM
, bGStat
, bNEMD
, bIterate
,
293 bFirstIterate
,bReadEkin
,bEkinAveVel
,bScaleEkin
, bConstrain
;
294 real ekin
,temp
,prescorr
,enercorr
,dvdlcorr
;
296 /* translate CGLO flags to booleans */
297 bRerunMD
= flags
& CGLO_RERUNMD
;
298 bStopCM
= flags
& CGLO_STOPCM
;
299 bGStat
= flags
& CGLO_GSTAT
;
300 bNEMD
= flags
& CGLO_NEMD
;
301 bReadEkin
= flags
& CGLO_READEKIN
;
302 bScaleEkin
= flags
& CGLO_SCALEEKIN
;
303 bEner
= flags
& CGLO_ENERGY
;
304 bTemp
= flags
& CGLO_TEMPERATURE
;
305 bPres
= flags
& CGLO_PRESSURE
;
306 bConstrain
= flags
& CGLO_CONSTRAINT
;
307 bIterate
= flags
& CGLO_ITERATE
;
308 bFirstIterate
= flags
& CGLO_FIRSTITERATE
;
310 /* we calculate a full state kinetic energy either with full-step velocity verlet
311 or half step where we need the pressure */
312 bEkinAveVel
= (ir
->eI
==eiVV
|| (ir
->eI
==eiVVAK
&& IR_NPT_TROTTER(ir
) && bPres
) || bReadEkin
);
314 /* in initalization, it sums the shake virial in vv, and to
315 sums ekinh_old in leapfrog (or if we are calculating ekinh_old for other reasons */
317 /* ########## Kinetic energy ############## */
321 /* Non-equilibrium MD: this is parallellized, but only does communication
322 * when there really is NEMD.
325 if (PAR(cr
) && (bNEMD
))
327 accumulate_u(cr
,&(ir
->opts
),ekind
);
332 restore_ekinstate_from_state(cr
,ekind
,&state_global
->ekinstate
);
337 calc_ke_part(state
,&(ir
->opts
),mdatoms
,ekind
,nrnb
,bEkinAveVel
,bIterate
);
342 /* Calculate center of mass velocity if necessary, also parallellized */
343 if (bStopCM
&& !bRerunMD
&& bEner
)
345 calc_vcm_grp(fplog
,mdatoms
->start
,mdatoms
->homenr
,mdatoms
,
346 state
->x
,state
->v
,vcm
);
350 if (bTemp
|| bPres
|| bEner
|| bConstrain
)
354 /* We will not sum ekinh_old,
355 * so signal that we still have to do it.
357 *bSumEkinhOld
= TRUE
;
364 for(i
=0; i
<eglsNR
; i
++)
366 gs_buf
[i
] = gs
->sig
[i
];
371 wallcycle_start(wcycle
,ewcMoveE
);
372 GMX_MPE_LOG(ev_global_stat_start
);
373 global_stat(fplog
,gstat
,cr
,enerd
,force_vir
,shake_vir
,mu_tot
,
375 gs
!= NULL
? eglsNR
: 0,gs_buf
,
377 *bSumEkinhOld
,flags
);
378 GMX_MPE_LOG(ev_global_stat_finish
);
379 wallcycle_stop(wcycle
,ewcMoveE
);
383 if (MULTISIM(cr
) && bInterSimGS
)
387 /* Communicate the signals between the simulations */
388 gmx_sum_sim(eglsNR
,gs_buf
,cr
->ms
);
390 /* Communicate the signals form the master to the others */
391 gmx_bcast(eglsNR
*sizeof(gs_buf
[0]),gs_buf
,cr
);
393 for(i
=0; i
<eglsNR
; i
++)
395 if (bInterSimGS
|| gs_simlocal
[i
])
397 /* Set the communicated signal only when it is non-zero,
398 * since signals might not be processed at each MD step.
400 gsi
= (gs_buf
[i
] >= 0 ?
401 (int)(gs_buf
[i
] + 0.5) :
402 (int)(gs_buf
[i
] - 0.5));
407 /* Turn off the local signal */
412 *bSumEkinhOld
= FALSE
;
416 if (!bNEMD
&& debug
&& bTemp
&& (vcm
->nr
> 0))
419 mdatoms
->start
,mdatoms
->start
+mdatoms
->homenr
,
420 state
->v
,vcm
->group_p
[0],
421 mdatoms
->massT
,mdatoms
->tmass
,ekind
->ekin
);
425 /* Do center of mass motion removal */
426 if (bStopCM
&& !bRerunMD
) /* is this correct? Does it get called too often with this logic? */
428 check_cm_grp(fplog
,vcm
,ir
,1);
429 do_stopcm_grp(fplog
,mdatoms
->start
,mdatoms
->homenr
,mdatoms
->cVCM
,
430 state
->x
,state
->v
,vcm
);
431 inc_nrnb(nrnb
,eNR_STOPCM
,mdatoms
->homenr
);
437 /* Sum the kinetic energies of the groups & calc temp */
438 /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with IR_NPT_TROTTER */
439 /* three maincase: VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
440 Leap with AveVel is also an option for the future but not supported now.
441 bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
442 If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
443 bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
444 If FALSE, we go ahead and erase over it.
446 enerd
->term
[F_TEMP
] = sum_ekin(&(ir
->opts
),ekind
,&(enerd
->term
[F_DKDL
]),
447 bEkinAveVel
,bIterate
,bScaleEkin
);
449 enerd
->term
[F_EKIN
] = trace(ekind
->ekin
);
452 /* ########## Long range energy information ###### */
454 if (bEner
|| bPres
|| bConstrain
)
456 calc_dispcorr(fplog
,ir
,fr
,0,top_global
->natoms
,box
,state
->lambda
,
457 corr_pres
,corr_vir
,&prescorr
,&enercorr
,&dvdlcorr
);
460 if (bEner
&& bFirstIterate
)
462 enerd
->term
[F_DISPCORR
] = enercorr
;
463 enerd
->term
[F_EPOT
] += enercorr
;
464 enerd
->term
[F_DVDL
] += dvdlcorr
;
465 if (fr
->efep
!= efepNO
) {
466 enerd
->dvdl_lin
+= dvdlcorr
;
470 /* ########## Now pressure ############## */
471 if (bPres
|| bConstrain
)
474 m_add(force_vir
,shake_vir
,total_vir
);
476 /* Calculate pressure and apply LR correction if PPPM is used.
477 * Use the box from last timestep since we already called update().
480 enerd
->term
[F_PRES
] = calc_pres(fr
->ePBC
,ir
->nwall
,box
,ekind
->ekin
,total_vir
,pres
,
481 (fr
->eeltype
==eelPPPM
)?enerd
->term
[F_COUL_RECIP
]:0.0);
483 /* Calculate long range corrections to pressure and energy */
484 /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
485 and computes enerd->term[F_DISPCORR]. Also modifies the
486 total_vir and pres tesors */
488 m_add(total_vir
,corr_vir
,total_vir
);
489 m_add(pres
,corr_pres
,pres
);
490 enerd
->term
[F_PDISPCORR
] = prescorr
;
491 enerd
->term
[F_PRES
] += prescorr
;
492 *pcurr
= enerd
->term
[F_PRES
];
493 /* calculate temperature using virial */
494 enerd
->term
[F_VTEMP
] = calc_temp(trace(total_vir
),ir
->opts
.nrdf
[0]);
500 /* Definitions for convergence of iterated constraints */
502 /* iterate constraints up to 50 times */
503 #define MAXITERCONST 50
508 real f
,fprev
,x
,xprev
;
511 real allrelerr
[MAXITERCONST
+2];
512 int num_close
; /* number of "close" violations, caused by limited precision. */
516 #define CONVERGEITER 0.000000001
517 #define CLOSE_ENOUGH 0.000001000
519 #define CONVERGEITER 0.0001
520 #define CLOSE_ENOUGH 0.0050
523 /* we want to keep track of the close calls. If there are too many, there might be some other issues.
524 so we make sure that it's either less than some predetermined number, or if more than that number,
525 only some small fraction of the total. */
526 #define MAX_NUMBER_CLOSE 50
527 #define FRACTION_CLOSE 0.001
529 /* maximum length of cyclic traps to check, emerging from limited numerical precision */
532 static void gmx_iterate_init(gmx_iterate_t
*iterate
,bool bIterate
)
537 iterate
->bIterate
= bIterate
;
538 iterate
->num_close
= 0;
539 for (i
=0;i
<MAXITERCONST
+2;i
++)
541 iterate
->allrelerr
[i
] = 0;
545 static bool done_iterating(const t_commrec
*cr
,FILE *fplog
, int nsteps
, gmx_iterate_t
*iterate
, bool bFirstIterate
, real fom
, real
*newf
)
547 /* monitor convergence, and use a secant search to propose new
550 The secant method computes x_{i+1} = x_{i} - f(x_{i}) * ---------------------
551 f(x_{i}) - f(x_{i-1})
553 The function we are trying to zero is fom-x, where fom is the
554 "figure of merit" which is the pressure (or the veta value) we
555 would get by putting in an old value of the pressure or veta into
556 the incrementor function for the step or half step. I have
557 verified that this gives the same answer as self consistent
558 iteration, usually in many fewer steps, especially for small tau_p.
560 We could possibly eliminate an iteration with proper use
561 of the value from the previous step, but that would take a bit
562 more bookkeeping, especially for veta, since tests indicate the
563 function of veta on the last step is not sufficiently close to
564 guarantee convergence this step. This is
565 good enough for now. On my tests, I could use tau_p down to
566 0.02, which is smaller that would ever be necessary in
567 practice. Generally, 3-5 iterations will be sufficient */
569 real relerr
,err
,xmin
;
577 iterate
->f
= fom
-iterate
->x
;
584 iterate
->f
= fom
-iterate
->x
; /* we want to zero this difference */
585 if ((iterate
->iter_i
> 1) && (iterate
->iter_i
< MAXITERCONST
))
587 if (iterate
->f
==iterate
->fprev
)
593 *newf
= iterate
->x
- (iterate
->x
-iterate
->xprev
)*(iterate
->f
)/(iterate
->f
-iterate
->fprev
);
598 /* just use self-consistent iteration the first step to initialize, or
599 if it's not converging (which happens occasionally -- need to investigate why) */
603 /* Consider a slight shortcut allowing us to exit one sooner -- we check the
604 difference between the closest of x and xprev to the new
605 value. To be 100% certain, we should check the difference between
606 the last result, and the previous result, or
608 relerr = (fabs((x-xprev)/fom));
610 but this is pretty much never necessary under typical conditions.
611 Checking numerically, it seems to lead to almost exactly the same
612 trajectories, but there are small differences out a few decimal
613 places in the pressure, and eventually in the v_eta, but it could
616 if (fabs(*newf-x) < fabs(*newf - xprev)) { xmin = x;} else { xmin = xprev;}
617 relerr = (fabs((*newf-xmin) / *newf));
620 err
= fabs((iterate
->f
-iterate
->fprev
));
621 relerr
= fabs(err
/fom
);
623 iterate
->allrelerr
[iterate
->iter_i
] = relerr
;
625 if (iterate
->iter_i
> 0)
629 fprintf(debug
,"Iterating NPT constraints: %6i %20.12f%14.6g%20.12f\n",
630 iterate
->iter_i
,fom
,relerr
,*newf
);
633 if ((relerr
< CONVERGEITER
) || (err
< CONVERGEITER
) || (fom
==0) || ((iterate
->x
== iterate
->xprev
) && iterate
->iter_i
> 1))
635 iterate
->bIterate
= FALSE
;
638 fprintf(debug
,"Iterating NPT constraints: CONVERGED\n");
642 if (iterate
->iter_i
> MAXITERCONST
)
644 if (relerr
< CLOSE_ENOUGH
)
647 for (i
=1;i
<CYCLEMAX
;i
++) {
648 if ((iterate
->allrelerr
[iterate
->iter_i
-(1+i
)] == iterate
->allrelerr
[iterate
->iter_i
-1]) &&
649 (iterate
->allrelerr
[iterate
->iter_i
-(1+i
)] == iterate
->allrelerr
[iterate
->iter_i
-(1+2*i
)])) {
653 fprintf(debug
,"Exiting from an NPT iterating cycle of length %d\n",i
);
660 /* step 1: trapped in a numerical attractor */
661 /* we are trapped in a numerical attractor, and can't converge any more, and are close to the final result.
662 Better to give up convergence here than have the simulation die.
664 iterate
->num_close
++;
669 /* Step #2: test if we are reasonably close for other reasons, then monitor the number. If not, die */
671 /* how many close calls have we had? If less than a few, we're OK */
672 if (iterate
->num_close
< MAX_NUMBER_CLOSE
)
674 sprintf(buf
,"Slight numerical convergence deviation with NPT at step %d, relative error only %10.5g, likely not a problem, continuing\n",nsteps
,relerr
);
675 md_print_warning(cr
,fplog
,buf
);
676 iterate
->num_close
++;
678 /* if more than a few, check the total fraction. If too high, die. */
679 } else if (iterate
->num_close
/(double)nsteps
> FRACTION_CLOSE
) {
680 gmx_fatal(FARGS
,"Could not converge NPT constraints, too many exceptions (%d%%\n",iterate
->num_close
/(double)nsteps
);
686 gmx_fatal(FARGS
,"Could not converge NPT constraints\n");
691 iterate
->xprev
= iterate
->x
;
693 iterate
->fprev
= iterate
->f
;
699 static void check_nst_param(FILE *fplog
,t_commrec
*cr
,
700 const char *desc_nst
,int nst
,
701 const char *desc_p
,int *p
)
705 if (*p
> 0 && *p
% nst
!= 0)
707 /* Round up to the next multiple of nst */
708 *p
= ((*p
)/nst
+ 1)*nst
;
709 sprintf(buf
,"NOTE: %s changes %s to %d\n",desc_nst
,desc_p
,*p
);
710 md_print_warning(cr
,fplog
,buf
);
714 static void reset_all_counters(FILE *fplog
,t_commrec
*cr
,
715 gmx_large_int_t step
,
716 gmx_large_int_t
*step_rel
,t_inputrec
*ir
,
717 gmx_wallcycle_t wcycle
,t_nrnb
*nrnb
,
718 gmx_runtime_t
*runtime
)
720 char buf
[STRLEN
],sbuf
[STEPSTRSIZE
];
722 /* Reset all the counters related to performance over the run */
723 sprintf(buf
,"Step %s: resetting all time and cycle counters\n",
724 gmx_step_str(step
,sbuf
));
725 md_print_warning(cr
,fplog
,buf
);
727 wallcycle_stop(wcycle
,ewcRUN
);
728 wallcycle_reset_all(wcycle
);
729 if (DOMAINDECOMP(cr
))
731 reset_dd_statistics_counters(cr
->dd
);
734 ir
->init_step
+= *step_rel
;
735 ir
->nsteps
-= *step_rel
;
737 wallcycle_start(wcycle
,ewcRUN
);
738 runtime_start(runtime
);
739 print_date_and_time(fplog
,cr
->nodeid
,"Restarted time",runtime
);
742 static int check_nstglobalcomm(FILE *fplog
,t_commrec
*cr
,
743 int nstglobalcomm
,t_inputrec
*ir
)
747 if (!EI_DYNAMICS(ir
->eI
))
752 if (nstglobalcomm
== -1)
754 if (ir
->nstcalcenergy
== 0 && ir
->nstlist
== 0)
757 if (ir
->nstenergy
> 0 && ir
->nstenergy
< nstglobalcomm
)
759 nstglobalcomm
= ir
->nstenergy
;
764 /* We assume that if nstcalcenergy > nstlist,
765 * nstcalcenergy is a multiple of nstlist.
767 if (ir
->nstcalcenergy
== 0 ||
768 (ir
->nstlist
> 0 && ir
->nstlist
< ir
->nstcalcenergy
))
770 nstglobalcomm
= ir
->nstlist
;
774 nstglobalcomm
= ir
->nstcalcenergy
;
780 if (ir
->nstlist
> 0 &&
781 nstglobalcomm
> ir
->nstlist
&& nstglobalcomm
% ir
->nstlist
!= 0)
783 nstglobalcomm
= (nstglobalcomm
/ ir
->nstlist
)*ir
->nstlist
;
784 sprintf(buf
,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm
);
785 md_print_warning(cr
,fplog
,buf
);
787 if (nstglobalcomm
> ir
->nstcalcenergy
)
789 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
790 "nstcalcenergy",&ir
->nstcalcenergy
);
793 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
794 "nstenergy",&ir
->nstenergy
);
796 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
797 "nstlog",&ir
->nstlog
);
800 if (ir
->comm_mode
!= ecmNO
&& ir
->nstcomm
< nstglobalcomm
)
802 sprintf(buf
,"WARNING: Changing nstcomm from %d to %d\n",
803 ir
->nstcomm
,nstglobalcomm
);
804 md_print_warning(cr
,fplog
,buf
);
805 ir
->nstcomm
= nstglobalcomm
;
808 return nstglobalcomm
;
811 void check_ir_old_tpx_versions(t_commrec
*cr
,FILE *fplog
,
812 t_inputrec
*ir
,gmx_mtop_t
*mtop
)
814 /* Check required for old tpx files */
815 if (IR_TWINRANGE(*ir
) && ir
->nstlist
> 1 &&
816 ir
->nstcalcenergy
% ir
->nstlist
!= 0)
818 md_print_warning(cr
,fplog
,"Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies");
820 if (gmx_mtop_ftype_count(mtop
,F_CONSTR
) +
821 gmx_mtop_ftype_count(mtop
,F_CONSTRNC
) > 0 &&
822 ir
->eConstrAlg
== econtSHAKE
)
824 md_print_warning(cr
,fplog
,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect");
825 if (ir
->epc
!= epcNO
)
827 gmx_fatal(FARGS
,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
830 check_nst_param(fplog
,cr
,"nstlist",ir
->nstlist
,
831 "nstcalcenergy",&ir
->nstcalcenergy
);
833 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
834 "nstenergy",&ir
->nstenergy
);
835 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
836 "nstlog",&ir
->nstlog
);
837 if (ir
->efep
!= efepNO
)
839 check_nst_param(fplog
,cr
,"nstdhdl",ir
->nstdhdl
,
840 "nstenergy",&ir
->nstenergy
);
845 bool bGStatEveryStep
;
846 gmx_large_int_t step_ns
;
847 gmx_large_int_t step_nscheck
;
858 static void reset_nlistheuristics(gmx_nlheur_t
*nlh
,gmx_large_int_t step
)
862 nlh
->step_nscheck
= step
;
865 static void init_nlistheuristics(gmx_nlheur_t
*nlh
,
866 bool bGStatEveryStep
,gmx_large_int_t step
)
868 nlh
->bGStatEveryStep
= bGStatEveryStep
;
875 reset_nlistheuristics(nlh
,step
);
878 static void update_nliststatistics(gmx_nlheur_t
*nlh
,gmx_large_int_t step
)
880 gmx_large_int_t nl_lt
;
881 char sbuf
[STEPSTRSIZE
],sbuf2
[STEPSTRSIZE
];
883 /* Determine the neighbor list life time */
884 nl_lt
= step
- nlh
->step_ns
;
887 fprintf(debug
,"%d atoms beyond ns buffer, updating neighbor list after %s steps\n",nlh
->nabnsb
,gmx_step_str(nl_lt
,sbuf
));
891 nlh
->s2
+= nl_lt
*nl_lt
;
892 nlh
->ab
+= nlh
->nabnsb
;
893 if (nlh
->lt_runav
== 0)
895 nlh
->lt_runav
= nl_lt
;
896 /* Initialize the fluctuation average
897 * such that at startup we check after 0 steps.
899 nlh
->lt_runav2
= sqr(nl_lt
/2.0);
901 /* Running average with 0.9 gives an exp. history of 9.5 */
902 nlh
->lt_runav2
= 0.9*nlh
->lt_runav2
+ 0.1*sqr(nlh
->lt_runav
- nl_lt
);
903 nlh
->lt_runav
= 0.9*nlh
->lt_runav
+ 0.1*nl_lt
;
904 if (nlh
->bGStatEveryStep
)
906 /* Always check the nlist validity */
907 nlh
->step_nscheck
= step
;
911 /* We check after: <life time> - 2*sigma
912 * The factor 2 is quite conservative,
913 * but we assume that with nstlist=-1 the user
914 * prefers exact integration over performance.
916 nlh
->step_nscheck
= step
917 + (int)(nlh
->lt_runav
- 2.0*sqrt(nlh
->lt_runav2
)) - 1;
921 fprintf(debug
,"nlist life time %s run av. %4.1f sig %3.1f check %s check with -gcom %d\n",
922 gmx_step_str(nl_lt
,sbuf
),nlh
->lt_runav
,sqrt(nlh
->lt_runav2
),
923 gmx_step_str(nlh
->step_nscheck
-step
+1,sbuf2
),
924 (int)(nlh
->lt_runav
- 2.0*sqrt(nlh
->lt_runav2
)));
928 static void set_nlistheuristics(gmx_nlheur_t
*nlh
,bool bReset
,gmx_large_int_t step
)
934 reset_nlistheuristics(nlh
,step
);
938 update_nliststatistics(nlh
,step
);
942 /* Initialize the cumulative coordinate scaling matrix */
943 clear_mat(nlh
->scale_tot
);
946 nlh
->scale_tot
[d
][d
] = 1.0;
950 double do_md(FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
951 const output_env_t oenv
, bool bVerbose
,bool bCompact
,
953 gmx_vsite_t
*vsite
,gmx_constr_t constr
,
954 int stepout
,t_inputrec
*ir
,
955 gmx_mtop_t
*top_global
,
957 t_state
*state_global
,
959 t_nrnb
*nrnb
,gmx_wallcycle_t wcycle
,
960 gmx_edsam_t ed
,t_forcerec
*fr
,
961 int repl_ex_nst
,int repl_ex_seed
,
962 real cpt_period
,real max_hours
,
963 const char *deviceOptions
,
965 gmx_runtime_t
*runtime
)
968 gmx_large_int_t step
,step_rel
;
971 bool bGStatEveryStep
,bGStat
,bNstEner
,bCalcEnerPres
;
972 bool bNS
,bNStList
,bSimAnn
,bStopCM
,bRerunMD
,bNotLastFrame
=FALSE
,
973 bFirstStep
,bStateFromTPX
,bInitStep
,bLastStep
,
974 bBornRadii
,bStartingFromCpt
;
976 bool bNEMD
,do_ene
,do_log
,do_verbose
,bRerunWarnNoV
=TRUE
,
977 bForceUpdate
=FALSE
,bCPT
;
980 int force_flags
,cglo_flags
;
981 tensor force_vir
,shake_vir
,total_vir
,tmp_vir
,pres
;
985 t_state
*bufstate
=NULL
;
986 matrix
*scale_tot
,pcoupl_mu
,M
,ebox
;
989 gmx_repl_ex_t repl_ex
=NULL
;
993 t_mdebin
*mdebin
=NULL
;
998 gmx_enerdata_t
*enerd
;
1000 gmx_global_stat_t gstat
;
1001 gmx_update_t upd
=NULL
;
1002 t_graph
*graph
=NULL
;
1006 gmx_groups_t
*groups
;
1007 gmx_ekindata_t
*ekind
, *ekind_save
;
1008 gmx_shellfc_t shellfc
;
1009 int count
,nconverged
=0;
1013 bool bTCR
=FALSE
,bConverged
=TRUE
,bOK
,bSumEkinhOld
,bExchanged
;
1015 bool bResetCountersHalfMaxH
=FALSE
;
1016 bool bVV
,bIterations
,bIterate
,bFirstIterate
,bTemp
,bPres
,bTrotter
;
1017 real temp0
,mu_aver
=0,dvdl
;
1019 atom_id
*grpindex
=NULL
;
1021 t_coupl_rec
*tcr
=NULL
;
1022 rvec
*xcopy
=NULL
,*vcopy
=NULL
,*cbuf
=NULL
;
1023 matrix boxcopy
={{0}},lastbox
;
1025 real fom
,oldfom
,veta_save
,pcurr
,scalevir
,tracevir
;
1028 real last_conserved
= 0;
1033 char sbuf
[STEPSTRSIZE
],sbuf2
[STEPSTRSIZE
];
1034 int handledSignal
=-1; /* compare to last_signal_recvd */
1035 gmx_iterate_t iterate
;
1037 /* Temporary addition for FAHCORE checkpointing */
1041 /* Check for special mdrun options */
1042 bRerunMD
= (Flags
& MD_RERUN
);
1043 bIonize
= (Flags
& MD_IONIZE
);
1044 bFFscan
= (Flags
& MD_FFSCAN
);
1045 bAppend
= (Flags
& MD_APPENDFILES
);
1046 if (Flags
& MD_RESETCOUNTERSHALFWAY
)
1050 /* Signal to reset the counters half the simulation steps. */
1051 wcycle_set_reset_counters(wcycle
,ir
->nsteps
/2);
1053 /* Signal to reset the counters halfway the simulation time. */
1054 bResetCountersHalfMaxH
= (max_hours
> 0);
1057 /* md-vv uses averaged full step velocities for T-control
1058 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
1059 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
1060 bVV
= EI_VV(ir
->eI
);
1061 if (bVV
) /* to store the initial velocities while computing virial */
1063 snew(cbuf
,top_global
->natoms
);
1065 /* all the iteratative cases - only if there are constraints */
1066 bIterations
= ((IR_NPT_TROTTER(ir
)) && (constr
) && (!bRerunMD
));
1067 bTrotter
= (bVV
&& (IR_NPT_TROTTER(ir
) || (IR_NVT_TROTTER(ir
))));
1071 /* Since we don't know if the frames read are related in any way,
1072 * rebuild the neighborlist at every step.
1075 ir
->nstcalcenergy
= 1;
1079 check_ir_old_tpx_versions(cr
,fplog
,ir
,top_global
);
1081 nstglobalcomm
= check_nstglobalcomm(fplog
,cr
,nstglobalcomm
,ir
);
1082 bGStatEveryStep
= (nstglobalcomm
== 1);
1084 if (!bGStatEveryStep
&& ir
->nstlist
== -1 && fplog
!= NULL
)
1087 "To reduce the energy communication with nstlist = -1\n"
1088 "the neighbor list validity should not be checked at every step,\n"
1089 "this means that exact integration is not guaranteed.\n"
1090 "The neighbor list validity is checked after:\n"
1091 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
1092 "In most cases this will result in exact integration.\n"
1093 "This reduces the energy communication by a factor of 2 to 3.\n"
1094 "If you want less energy communication, set nstlist > 3.\n\n");
1097 if (bRerunMD
|| bFFscan
)
1101 groups
= &top_global
->groups
;
1103 /* Initial values */
1104 init_md(fplog
,cr
,ir
,oenv
,&t
,&t0
,&state_global
->lambda
,&lam0
,
1105 nrnb
,top_global
,&upd
,
1106 nfile
,fnm
,&outf
,&mdebin
,
1107 force_vir
,shake_vir
,mu_tot
,&bNEMD
,&bSimAnn
,&vcm
,state_global
,Flags
);
1109 clear_mat(total_vir
);
1111 /* Energy terms and groups */
1113 init_enerdata(top_global
->groups
.grps
[egcENER
].nr
,ir
->n_flambda
,enerd
);
1114 if (DOMAINDECOMP(cr
))
1120 snew(f
,top_global
->natoms
);
1123 /* Kinetic energy data */
1125 init_ekindata(fplog
,top_global
,&(ir
->opts
),ekind
);
1126 /* needed for iteration of constraints */
1128 init_ekindata(fplog
,top_global
,&(ir
->opts
),ekind_save
);
1129 /* Copy the cos acceleration to the groups struct */
1130 ekind
->cosacc
.cos_accel
= ir
->cos_accel
;
1132 gstat
= global_stat_init(ir
);
1135 /* Check for polarizable models and flexible constraints */
1136 shellfc
= init_shell_flexcon(fplog
,
1137 top_global
,n_flexible_constraints(constr
),
1138 (ir
->bContinuation
||
1139 (DOMAINDECOMP(cr
) && !MASTER(cr
))) ?
1140 NULL
: state_global
->x
);
1145 tMPI_Thread_mutex_lock(&deform_init_box_mutex
);
1147 set_deform_reference_box(upd
,
1148 deform_init_init_step_tpx
,
1149 deform_init_box_tpx
);
1151 tMPI_Thread_mutex_unlock(&deform_init_box_mutex
);
1156 double io
= compute_io(ir
,top_global
->natoms
,groups
,mdebin
->ebin
->nener
,1);
1157 if ((io
> 2000) && MASTER(cr
))
1159 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
1163 if (DOMAINDECOMP(cr
)) {
1164 top
= dd_init_local_top(top_global
);
1167 dd_init_local_state(cr
->dd
,state_global
,state
);
1169 if (DDMASTER(cr
->dd
) && ir
->nstfout
) {
1170 snew(f_global
,state_global
->natoms
);
1174 /* Initialize the particle decomposition and split the topology */
1175 top
= split_system(fplog
,top_global
,ir
,cr
);
1177 pd_cg_range(cr
,&fr
->cg0
,&fr
->hcg
);
1178 pd_at_range(cr
,&a0
,&a1
);
1180 top
= gmx_mtop_generate_local_top(top_global
,ir
);
1183 a1
= top_global
->natoms
;
1186 state
= partdec_init_local_state(cr
,state_global
);
1189 atoms2md(top_global
,ir
,0,NULL
,a0
,a1
-a0
,mdatoms
);
1192 set_vsite_top(vsite
,top
,mdatoms
,cr
);
1195 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
) {
1196 graph
= mk_graph(fplog
,&(top
->idef
),0,top_global
->natoms
,FALSE
,FALSE
);
1200 make_local_shells(cr
,mdatoms
,shellfc
);
1203 if (ir
->pull
&& PAR(cr
)) {
1204 dd_make_local_pull_groups(NULL
,ir
->pull
,mdatoms
);
1208 if (DOMAINDECOMP(cr
))
1210 /* Distribute the charge groups over the nodes from the master node */
1211 dd_partition_system(fplog
,ir
->init_step
,cr
,TRUE
,1,
1212 state_global
,top_global
,ir
,
1213 state
,&f
,mdatoms
,top
,fr
,
1214 vsite
,shellfc
,constr
,
1218 update_mdatoms(mdatoms
,state
->lambda
);
1222 /* Update mdebin with energy history if appending to output files */
1223 if ( Flags
& MD_APPENDFILES
)
1225 restore_energyhistory_from_state(mdebin
,&state_global
->enerhist
);
1227 /* Set the initial energy history in state to zero by updating once */
1228 update_energyhistory(&state_global
->enerhist
,mdebin
);
1231 if ((state
->flags
& (1<<estLD_RNG
)) && (Flags
& MD_READ_RNG
)) {
1232 /* Set the random state if we read a checkpoint file */
1233 set_stochd_state(upd
,state
);
1236 /* Initialize constraints */
1238 if (!DOMAINDECOMP(cr
))
1239 set_constraints(constr
,top
,ir
,mdatoms
,cr
);
1242 /* Check whether we have to GCT stuff */
1243 bTCR
= ftp2bSet(efGCT
,nfile
,fnm
);
1246 fprintf(stderr
,"Will do General Coupling Theory!\n");
1248 gnx
= top_global
->mols
.nr
;
1250 for(i
=0; (i
<gnx
); i
++) {
1255 if (repl_ex_nst
> 0 && MASTER(cr
))
1256 repl_ex
= init_replica_exchange(fplog
,cr
->ms
,state_global
,ir
,
1257 repl_ex_nst
,repl_ex_seed
);
1259 if (!ir
->bContinuation
&& !bRerunMD
)
1261 if (mdatoms
->cFREEZE
&& (state
->flags
& (1<<estV
)))
1263 /* Set the velocities of frozen particles to zero */
1264 for(i
=mdatoms
->start
; i
<mdatoms
->start
+mdatoms
->homenr
; i
++)
1266 for(m
=0; m
<DIM
; m
++)
1268 if (ir
->opts
.nFreeze
[mdatoms
->cFREEZE
[i
]][m
])
1278 /* Constrain the initial coordinates and velocities */
1279 do_constrain_first(fplog
,constr
,ir
,mdatoms
,state
,f
,
1280 graph
,cr
,nrnb
,fr
,top
,shake_vir
);
1284 /* Construct the virtual sites for the initial configuration */
1285 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,NULL
,
1286 top
->idef
.iparams
,top
->idef
.il
,
1287 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
1293 /* I'm assuming we need global communication the first time! MRS */
1294 cglo_flags
= (CGLO_TEMPERATURE
| CGLO_GSTAT
1295 | (bVV
? CGLO_PRESSURE
:0)
1296 | (bVV
? CGLO_CONSTRAINT
:0)
1297 | (bRerunMD
? CGLO_RERUNMD
:0)
1298 | ((Flags
& MD_READ_EKIN
) ? CGLO_READEKIN
:0));
1300 bSumEkinhOld
= FALSE
;
1301 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
1302 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
1303 constr
,NULL
,FALSE
,state
->box
,
1304 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,cglo_flags
);
1305 if (ir
->eI
== eiVVAK
) {
1306 /* a second call to get the half step temperature initialized as well */
1307 /* we do the same call as above, but turn the pressure off -- internally, this
1308 is recognized as a velocity verlet half-step kinetic energy calculation.
1309 This minimized excess variables, but perhaps loses some logic?*/
1311 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
1312 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
1313 constr
,NULL
,FALSE
,state
->box
,
1314 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
1315 cglo_flags
&~ CGLO_PRESSURE
);
1318 /* Calculate the initial half step temperature, and save the ekinh_old */
1319 if (!(Flags
& MD_STARTFROMCPT
))
1321 for(i
=0; (i
<ir
->opts
.ngtc
); i
++)
1323 copy_mat(ekind
->tcstat
[i
].ekinh
,ekind
->tcstat
[i
].ekinh_old
);
1326 temp0
= enerd
->term
[F_TEMP
];
1328 /* if using an iterative algorithm, we need to create a working directory for the state. */
1331 bufstate
= init_bufstate(state
);
1335 snew(xcopy
,state
->natoms
);
1336 snew(vcopy
,state
->natoms
);
1337 copy_rvecn(state
->x
,xcopy
,0,state
->natoms
);
1338 copy_rvecn(state
->v
,vcopy
,0,state
->natoms
);
1339 copy_mat(state
->box
,boxcopy
);
1342 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
1343 temperature control */
1344 trotter_seq
= init_npt_vars(ir
,state
,&MassQ
,bTrotter
);
1348 if (constr
&& !ir
->bContinuation
&& ir
->eConstrAlg
== econtLINCS
)
1351 "RMS relative constraint deviation after constraining: %.2e\n",
1352 constr_rmsd(constr
,FALSE
));
1356 enerd
->term
[F_TEMP
] *= 2; /* result of averages being done over previous and current step,
1357 and there is no previous step */
1359 fprintf(fplog
,"Initial temperature: %g K\n",enerd
->term
[F_TEMP
]);
1362 fprintf(stderr
,"starting md rerun '%s', reading coordinates from"
1363 " input trajectory '%s'\n\n",
1364 *(top_global
->name
),opt2fn("-rerun",nfile
,fnm
));
1367 fprintf(stderr
,"Calculated time to finish depends on nsteps from "
1368 "run input file,\nwhich may not correspond to the time "
1369 "needed to process input trajectory.\n\n");
1375 fprintf(stderr
,"starting mdrun '%s'\n",
1376 *(top_global
->name
));
1377 if (ir
->nsteps
>= 0)
1379 sprintf(tbuf
,"%8.1f",(ir
->init_step
+ir
->nsteps
)*ir
->delta_t
);
1383 sprintf(tbuf
,"%s","infinite");
1385 if (ir
->init_step
> 0)
1387 fprintf(stderr
,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
1388 gmx_step_str(ir
->init_step
+ir
->nsteps
,sbuf
),tbuf
,
1389 gmx_step_str(ir
->init_step
,sbuf2
),
1390 ir
->init_step
*ir
->delta_t
);
1394 fprintf(stderr
,"%s steps, %s ps.\n",
1395 gmx_step_str(ir
->nsteps
,sbuf
),tbuf
);
1398 fprintf(fplog
,"\n");
1401 /* Set and write start time */
1402 runtime_start(runtime
);
1403 print_date_and_time(fplog
,cr
->nodeid
,"Started mdrun",runtime
);
1404 wallcycle_start(wcycle
,ewcRUN
);
1406 fprintf(fplog
,"\n");
1408 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
1410 chkpt_ret
=fcCheckPointParallel( cr
->nodeid
,
1412 if ( chkpt_ret
== 0 )
1413 gmx_fatal( 3,__FILE__
,__LINE__
, "Checkpoint error on step %d\n", 0 );
1417 /***********************************************************
1419 * Loop over MD steps
1421 ************************************************************/
1423 /* if rerunMD then read coordinates and velocities from input trajectory */
1426 if (getenv("GMX_FORCE_UPDATE"))
1428 bForceUpdate
= TRUE
;
1431 bNotLastFrame
= read_first_frame(oenv
,&status
,
1432 opt2fn("-rerun",nfile
,fnm
),
1433 &rerun_fr
,TRX_NEED_X
| TRX_READ_V
);
1434 if (rerun_fr
.natoms
!= top_global
->natoms
)
1437 "Number of atoms in trajectory (%d) does not match the "
1438 "run input file (%d)\n",
1439 rerun_fr
.natoms
,top_global
->natoms
);
1441 if (ir
->ePBC
!= epbcNONE
)
1445 gmx_fatal(FARGS
,"Rerun trajectory frame step %d time %f does not contain a box, while pbc is used",rerun_fr
.step
,rerun_fr
.time
);
1447 if (max_cutoff2(ir
->ePBC
,rerun_fr
.box
) < sqr(fr
->rlistlong
))
1449 gmx_fatal(FARGS
,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr
.step
,rerun_fr
.time
);
1452 /* Set the shift vectors.
1453 * Necessary here when have a static box different from the tpr box.
1455 calc_shifts(rerun_fr
.box
,fr
->shift_vec
);
1459 /* loop over MD steps or if rerunMD to end of input trajectory */
1461 /* Skip the first Nose-Hoover integration when we get the state from tpx */
1462 bStateFromTPX
= !opt2bSet("-cpi",nfile
,fnm
);
1463 bInitStep
= bFirstStep
&& (bStateFromTPX
|| bVV
);
1464 bStartingFromCpt
= (Flags
& MD_STARTFROMCPT
) && bInitStep
;
1466 bSumEkinhOld
= FALSE
;
1469 init_global_signals(&gs
,cr
,ir
,repl_ex_nst
);
1471 step
= ir
->init_step
;
1474 if (ir
->nstlist
== -1)
1476 init_nlistheuristics(&nlh
,bGStatEveryStep
,step
);
1479 bLastStep
= (bRerunMD
|| (ir
->nsteps
>= 0 && step_rel
> ir
->nsteps
));
1480 while (!bLastStep
|| (bRerunMD
&& bNotLastFrame
)) {
1482 wallcycle_start(wcycle
,ewcSTEP
);
1484 GMX_MPE_LOG(ev_timestep1
);
1487 if (rerun_fr
.bStep
) {
1488 step
= rerun_fr
.step
;
1489 step_rel
= step
- ir
->init_step
;
1491 if (rerun_fr
.bTime
) {
1501 bLastStep
= (step_rel
== ir
->nsteps
);
1502 t
= t0
+ step
*ir
->delta_t
;
1505 if (ir
->efep
!= efepNO
)
1507 if (bRerunMD
&& rerun_fr
.bLambda
&& (ir
->delta_lambda
!=0))
1509 state_global
->lambda
= rerun_fr
.lambda
;
1513 state_global
->lambda
= lam0
+ step
*ir
->delta_lambda
;
1515 state
->lambda
= state_global
->lambda
;
1516 bDoDHDL
= do_per_step(step
,ir
->nstdhdl
);
1521 update_annealing_target_temp(&(ir
->opts
),t
);
1526 if (!(DOMAINDECOMP(cr
) && !MASTER(cr
)))
1528 for(i
=0; i
<state_global
->natoms
; i
++)
1530 copy_rvec(rerun_fr
.x
[i
],state_global
->x
[i
]);
1534 for(i
=0; i
<state_global
->natoms
; i
++)
1536 copy_rvec(rerun_fr
.v
[i
],state_global
->v
[i
]);
1541 for(i
=0; i
<state_global
->natoms
; i
++)
1543 clear_rvec(state_global
->v
[i
]);
1547 fprintf(stderr
,"\nWARNING: Some frames do not contain velocities.\n"
1548 " Ekin, temperature and pressure are incorrect,\n"
1549 " the virial will be incorrect when constraints are present.\n"
1551 bRerunWarnNoV
= FALSE
;
1555 copy_mat(rerun_fr
.box
,state_global
->box
);
1556 copy_mat(state_global
->box
,state
->box
);
1558 if (vsite
&& (Flags
& MD_RERUN_VSITE
))
1560 if (DOMAINDECOMP(cr
))
1562 gmx_fatal(FARGS
,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
1566 /* Following is necessary because the graph may get out of sync
1567 * with the coordinates if we only have every N'th coordinate set
1569 mk_mshift(fplog
,graph
,fr
->ePBC
,state
->box
,state
->x
);
1570 shift_self(graph
,state
->box
,state
->x
);
1572 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,state
->v
,
1573 top
->idef
.iparams
,top
->idef
.il
,
1574 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
1577 unshift_self(graph
,state
->box
,state
->x
);
1582 /* Stop Center of Mass motion */
1583 bStopCM
= (ir
->comm_mode
!= ecmNO
&& do_per_step(step
,ir
->nstcomm
));
1585 /* Copy back starting coordinates in case we're doing a forcefield scan */
1588 for(ii
=0; (ii
<state
->natoms
); ii
++)
1590 copy_rvec(xcopy
[ii
],state
->x
[ii
]);
1591 copy_rvec(vcopy
[ii
],state
->v
[ii
]);
1593 copy_mat(boxcopy
,state
->box
);
1598 /* for rerun MD always do Neighbour Searching */
1599 bNS
= (bFirstStep
|| ir
->nstlist
!= 0);
1604 /* Determine whether or not to do Neighbour Searching and LR */
1605 bNStList
= (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0);
1607 bNS
= (bFirstStep
|| bExchanged
|| bNStList
||
1608 (ir
->nstlist
== -1 && nlh
.nabnsb
> 0));
1610 if (bNS
&& ir
->nstlist
== -1)
1612 set_nlistheuristics(&nlh
,bFirstStep
|| bExchanged
,step
);
1616 if (gs
.set
[eglsTERM
] > 0 || (gs
.set
[eglsTERM
] < 0 && bNS
))
1621 /* Determine whether or not to update the Born radii if doing GB */
1622 bBornRadii
=bFirstStep
;
1623 if (ir
->implicit_solvent
&& (step
% ir
->nstgbradii
==0))
1628 do_log
= do_per_step(step
,ir
->nstlog
) || bFirstStep
|| bLastStep
;
1629 do_verbose
= bVerbose
&&
1630 (step
% stepout
== 0 || bFirstStep
|| bLastStep
);
1632 if (bNS
&& !(bFirstStep
&& ir
->bContinuation
&& !bRerunMD
))
1636 bMasterState
= TRUE
;
1640 bMasterState
= FALSE
;
1641 /* Correct the new box if it is too skewed */
1642 if (DYNAMIC_BOX(*ir
))
1644 if (correct_box(fplog
,step
,state
->box
,graph
))
1646 bMasterState
= TRUE
;
1649 if (DOMAINDECOMP(cr
) && bMasterState
)
1651 dd_collect_state(cr
->dd
,state
,state_global
);
1655 if (DOMAINDECOMP(cr
))
1657 /* Repartition the domain decomposition */
1658 wallcycle_start(wcycle
,ewcDOMDEC
);
1659 dd_partition_system(fplog
,step
,cr
,
1660 bMasterState
,nstglobalcomm
,
1661 state_global
,top_global
,ir
,
1662 state
,&f
,mdatoms
,top
,fr
,
1663 vsite
,shellfc
,constr
,
1664 nrnb
,wcycle
,do_verbose
);
1665 wallcycle_stop(wcycle
,ewcDOMDEC
);
1666 /* If using an iterative integrator, reallocate space to match the decomposition */
1670 if (MASTER(cr
) && do_log
&& !bFFscan
)
1672 print_ebin_header(fplog
,step
,t
,state
->lambda
);
1675 if (ir
->efep
!= efepNO
)
1677 update_mdatoms(mdatoms
,state
->lambda
);
1680 if (bRerunMD
&& rerun_fr
.bV
)
1683 /* We need the kinetic energy at minus the half step for determining
1684 * the full step kinetic energy and possibly for T-coupling.*/
1685 /* This may not be quite working correctly yet . . . . */
1686 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
1687 wcycle
,enerd
,NULL
,NULL
,NULL
,NULL
,mu_tot
,
1688 constr
,NULL
,FALSE
,state
->box
,
1689 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
1690 CGLO_RERUNMD
| CGLO_GSTAT
| CGLO_TEMPERATURE
);
1692 clear_mat(force_vir
);
1694 /* Ionize the atoms if necessary */
1697 ionize(fplog
,oenv
,mdatoms
,top_global
,t
,ir
,state
->x
,state
->v
,
1698 mdatoms
->start
,mdatoms
->start
+mdatoms
->homenr
,state
->box
,cr
);
1701 /* Update force field in ffscan program */
1704 if (update_forcefield(fplog
,
1706 mdatoms
->nr
,state
->x
,state
->box
)) {
1707 if (gmx_parallel_env_initialized())
1715 GMX_MPE_LOG(ev_timestep2
);
1717 /* We write a checkpoint at this MD step when:
1718 * either at an NS step when we signalled through gs,
1719 * or at the last step (but not when we do not want confout),
1720 * but never at the first step or with rerun.
1722 bCPT
= (((gs
.set
[eglsCHKPT
] && bNS
) ||
1723 (bLastStep
&& (Flags
& MD_CONFOUT
))) &&
1724 step
> ir
->init_step
&& !bRerunMD
);
1727 gs
.set
[eglsCHKPT
] = 0;
1730 /* Determine the energy and pressure:
1731 * at nstcalcenergy steps and at energy output steps (set below).
1733 bNstEner
= (bGStatEveryStep
|| do_per_step(step
,ir
->nstcalcenergy
));
1734 bCalcEnerPres
= bNstEner
;
1736 /* Do we need global communication ? */
1737 bGStat
= (bCalcEnerPres
|| bStopCM
||
1738 (ir
->nstlist
== -1 && !bRerunMD
&& step
>= nlh
.step_nscheck
));
1740 do_ene
= (do_per_step(step
,ir
->nstenergy
) || bLastStep
);
1742 if (do_ene
|| do_log
)
1744 bCalcEnerPres
= TRUE
;
1748 /* these CGLO_ options remain the same throughout the iteration */
1749 cglo_flags
= ((bRerunMD
? CGLO_RERUNMD
: 0) |
1750 (bStopCM
? CGLO_STOPCM
: 0) |
1751 (bGStat
? CGLO_GSTAT
: 0) |
1752 (bNEMD
? CGLO_NEMD
: 0)
1755 force_flags
= (GMX_FORCE_STATECHANGED
|
1756 ((DYNAMIC_BOX(*ir
) || bRerunMD
) ? GMX_FORCE_DYNAMICBOX
: 0) |
1757 GMX_FORCE_ALLFORCES
|
1758 (bNStList
? GMX_FORCE_DOLR
: 0) |
1760 (bCalcEnerPres
? GMX_FORCE_VIRIAL
: 0) |
1761 (bDoDHDL
? GMX_FORCE_DHDL
: 0)
1766 /* Now is the time to relax the shells */
1767 count
=relax_shell_flexcon(fplog
,cr
,bVerbose
,bFFscan
? step
+1 : step
,
1769 bStopCM
,top
,top_global
,
1771 state
,f
,force_vir
,mdatoms
,
1772 nrnb
,wcycle
,graph
,groups
,
1773 shellfc
,fr
,bBornRadii
,t
,mu_tot
,
1774 state
->natoms
,&bConverged
,vsite
,
1785 /* The coordinates (x) are shifted (to get whole molecules)
1787 * This is parallellized as well, and does communication too.
1788 * Check comments in sim_util.c
1791 do_force(fplog
,cr
,ir
,step
,nrnb
,wcycle
,top
,top_global
,groups
,
1792 state
->box
,state
->x
,&state
->hist
,
1793 f
,force_vir
,mdatoms
,enerd
,fcd
,
1794 state
->lambda
,graph
,
1795 fr
,vsite
,mu_tot
,t
,outf
->fp_field
,ed
,bBornRadii
,
1796 (bNS
? GMX_FORCE_NS
: 0) | force_flags
);
1799 GMX_BARRIER(cr
->mpi_comm_mygroup
);
1803 mu_aver
= calc_mu_aver(cr
,state
->x
,mdatoms
->chargeA
,
1804 mu_tot
,&top_global
->mols
,mdatoms
,gnx
,grpindex
);
1807 if (bTCR
&& bFirstStep
)
1809 tcr
=init_coupling(fplog
,nfile
,fnm
,cr
,fr
,mdatoms
,&(top
->idef
));
1810 fprintf(fplog
,"Done init_coupling\n");
1814 /* ############### START FIRST UPDATE HALF-STEP ############### */
1816 if (!bStartingFromCpt
&& !bRerunMD
)
1822 /* if using velocity verlet with full time step Ekin,
1823 * take the first half step only to compute the
1824 * virial for the first step. From there,
1825 * revert back to the initial coordinates
1826 * so that the input is actually the initial step.
1828 copy_rvecn(state
->v
,cbuf
,0,state
->natoms
); /* should make this better for parallelizing? */
1831 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1834 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[1]);
1837 update_coords(fplog
,step
,ir
,mdatoms
,state
,
1838 f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
1839 ekind
,M
,wcycle
,upd
,bInitStep
,etrtVELOCITY
,
1840 cr
,nrnb
,constr
,&top
->idef
);
1844 gmx_iterate_init(&iterate
,bIterations
&& !bInitStep
);
1846 /* for iterations, we save these vectors, as we will be self-consistently iterating
1848 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1850 /* save the state */
1851 if (bIterations
&& iterate
.bIterate
) {
1852 copy_coupling_state(state
,bufstate
,ekind
,ekind_save
,&(ir
->opts
));
1856 bFirstIterate
= TRUE
;
1857 while (bFirstIterate
|| (bIterations
&& iterate
.bIterate
))
1859 if (bIterations
&& iterate
.bIterate
)
1861 copy_coupling_state(bufstate
,state
,ekind_save
,ekind
,&(ir
->opts
));
1862 if (bFirstIterate
&& bTrotter
)
1864 /* The first time through, we need a decent first estimate
1865 of veta(t+dt) to compute the constraints. Do
1866 this by computing the box volume part of the
1867 trotter integration at this time. Nothing else
1868 should be changed by this routine here. If
1869 !(first time), we start with the previous value
1872 veta_save
= state
->veta
;
1873 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[0]);
1874 vetanew
= state
->veta
;
1875 state
->veta
= veta_save
;
1880 if ( !bRerunMD
|| rerun_fr
.bV
|| bForceUpdate
) { /* Why is rerun_fr.bV here? Unclear. */
1883 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
1884 &top
->idef
,shake_vir
,NULL
,
1885 cr
,nrnb
,wcycle
,upd
,constr
,
1886 bInitStep
,TRUE
,bCalcEnerPres
,vetanew
);
1888 if (!bOK
&& !bFFscan
)
1890 gmx_fatal(FARGS
,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1895 { /* Need to unshift here if a do_force has been
1896 called in the previous step */
1897 unshift_self(graph
,state
->box
,state
->x
);
1902 /* if VV, compute the pressure and constraints */
1903 /* if VV2, the pressure and constraints only if using pressure control.*/
1904 bPres
= (ir
->eI
==eiVV
|| IR_NPT_TROTTER(ir
));
1905 bTemp
= ((ir
->eI
==eiVV
&&(!bInitStep
)) || (ir
->eI
==eiVVAK
&& IR_NPT_TROTTER(ir
)));
1906 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
1907 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
1908 constr
,NULL
,FALSE
,state
->box
,
1909 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
1912 | (bTemp
? CGLO_TEMPERATURE
:0)
1913 | (bPres
? CGLO_PRESSURE
: 0)
1914 | (bPres
? CGLO_CONSTRAINT
: 0)
1915 | (iterate
.bIterate
? CGLO_ITERATE
: 0)
1916 | (bFirstIterate
? CGLO_FIRSTITERATE
: 0)
1920 /* explanation of above:
1921 a) We compute Ekin at the full time step
1922 if 1) we are using the AveVel Ekin, and it's not the
1923 initial step, or 2) if we are using AveEkin, but need the full
1924 time step kinetic energy for the pressure.
1925 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
1926 EkinAveVel because it's needed for the pressure */
1928 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1929 if (bVV
&& !bInitStep
)
1931 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[2]);
1935 done_iterating(cr
,fplog
,step
,&iterate
,bFirstIterate
,
1936 state
->veta
,&vetanew
))
1940 bFirstIterate
= FALSE
;
1943 if (bTrotter
&& !bInitStep
) {
1944 copy_mat(shake_vir
,state
->svir_prev
);
1945 copy_mat(force_vir
,state
->fvir_prev
);
1946 if (IR_NVT_TROTTER(ir
) && ir
->eI
==eiVV
) {
1947 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1948 enerd
->term
[F_TEMP
] = sum_ekin(&(ir
->opts
),ekind
,NULL
,(ir
->eI
==eiVV
),FALSE
,FALSE
);
1949 enerd
->term
[F_EKIN
] = trace(ekind
->ekin
);
1952 /* if it's the initial step, we performed this first step just to get the constraint virial */
1953 if (bInitStep
&& ir
->eI
==eiVV
) {
1954 copy_rvecn(cbuf
,state
->v
,0,state
->natoms
);
1957 if (fr
->bSepDVDL
&& fplog
&& do_log
)
1959 fprintf(fplog
,sepdvdlformat
,"Constraint",0.0,dvdl
);
1961 enerd
->term
[F_DHDL_CON
] += dvdl
;
1963 GMX_MPE_LOG(ev_timestep1
);
1967 /* MRS -- now done iterating -- compute the conserved quantity */
1970 if (IR_NVT_TROTTER(ir
) || IR_NPT_TROTTER(ir
))
1973 NPT_energy(ir
,state
,&MassQ
);
1974 if ((ir
->eDispCorr
!= edispcEnerPres
) && (ir
->eDispCorr
!= edispcAllEnerPres
))
1976 last_conserved
-= enerd
->term
[F_DISPCORR
];
1980 last_ekin
= enerd
->term
[F_EKIN
]; /* does this get preserved through checkpointing? */
1984 /* ######## END FIRST UPDATE STEP ############## */
1985 /* ######## If doing VV, we now have v(dt) ###### */
1987 /* ################## START TRAJECTORY OUTPUT ################# */
1989 /* Now we have the energies and forces corresponding to the
1990 * coordinates at time t. We must output all of this before
1992 * for RerunMD t is read from input trajectory
1994 GMX_MPE_LOG(ev_output_start
);
1997 if (do_per_step(step
,ir
->nstxout
)) { mdof_flags
|= MDOF_X
; }
1998 if (do_per_step(step
,ir
->nstvout
)) { mdof_flags
|= MDOF_V
; }
1999 if (do_per_step(step
,ir
->nstfout
)) { mdof_flags
|= MDOF_F
; }
2000 if (do_per_step(step
,ir
->nstxtcout
)) { mdof_flags
|= MDOF_XTC
; }
2001 if (bCPT
) { mdof_flags
|= MDOF_CPT
; };
2005 fcReportProgress( ir
->nsteps
, step
);
2009 /* Enforce writing positions and velocities at end of run */
2010 mdof_flags
|= (MDOF_X
| MDOF_V
);
2013 int nthreads
=(cr
->nthreads
==0 ? 1 : cr
->nthreads
);
2014 int nnodes
=(cr
->nnodes
==0 ? 1 : cr
->nnodes
);
2016 /*Gromacs drives checkpointing; no ||
2017 fcCheckPointPendingThreads(cr->nodeid,
2019 /* sync bCPT and fc record-keeping */
2020 if (bCPT
&& MASTER(cr
))
2021 fcRequestCheckPoint();
2025 if (mdof_flags
!= 0)
2027 wallcycle_start(wcycle
,ewcTRAJ
);
2030 if (state
->flags
& (1<<estLD_RNG
))
2032 get_stochd_state(upd
,state
);
2038 state_global
->ekinstate
.bUpToDate
= FALSE
;
2042 update_ekinstate(&state_global
->ekinstate
,ekind
);
2043 state_global
->ekinstate
.bUpToDate
= TRUE
;
2045 update_energyhistory(&state_global
->enerhist
,mdebin
);
2048 write_traj(fplog
,cr
,outf
,mdof_flags
,top_global
,
2049 step
,t
,state
,state_global
,f
,f_global
,&n_xtc
,&x_xtc
);
2056 if (bLastStep
&& step_rel
== ir
->nsteps
&&
2057 (Flags
& MD_CONFOUT
) && MASTER(cr
) &&
2058 !bRerunMD
&& !bFFscan
)
2060 /* x and v have been collected in write_traj,
2061 * because a checkpoint file will always be written
2064 fprintf(stderr
,"\nWriting final coordinates.\n");
2065 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
&&
2068 /* Make molecules whole only for confout writing */
2069 do_pbc_mtop(fplog
,ir
->ePBC
,state
->box
,top_global
,state_global
->x
);
2071 write_sto_conf_mtop(ftp2fn(efSTO
,nfile
,fnm
),
2072 *top_global
->name
,top_global
,
2073 state_global
->x
,state_global
->v
,
2074 ir
->ePBC
,state
->box
);
2077 wallcycle_stop(wcycle
,ewcTRAJ
);
2079 GMX_MPE_LOG(ev_output_finish
);
2081 /* kludge -- virial is lost with restart for NPT control. Must restart */
2082 if (bStartingFromCpt
&& bVV
)
2084 copy_mat(state
->svir_prev
,shake_vir
);
2085 copy_mat(state
->fvir_prev
,force_vir
);
2087 /* ################## END TRAJECTORY OUTPUT ################ */
2089 /* Determine the pressure:
2090 * always when we want exact averages in the energy file,
2091 * at ns steps when we have pressure coupling,
2092 * otherwise only at energy output steps (set below).
2095 bNstEner
= (bGStatEveryStep
|| do_per_step(step
,ir
->nstcalcenergy
));
2096 bCalcEnerPres
= bNstEner
;
2098 /* Do we need global communication ? */
2099 bGStat
= (bGStatEveryStep
|| bStopCM
|| bNS
||
2100 (ir
->nstlist
== -1 && !bRerunMD
&& step
>= nlh
.step_nscheck
));
2102 do_ene
= (do_per_step(step
,ir
->nstenergy
) || bLastStep
);
2104 if (do_ene
|| do_log
)
2106 bCalcEnerPres
= TRUE
;
2110 /* Determine the wallclock run time up till now */
2111 run_time
= gmx_gettime() - (double)runtime
->real
;
2113 /* Check whether everything is still allright */
2114 if ((bGotStopNextStepSignal
|| bGotStopNextNSStepSignal
) &&
2115 (handledSignal
!=last_signal_number_recvd
) &&
2118 if (bGotStopNextStepSignal
|| ir
->nstlist
== 0)
2120 gs
.sig
[eglsTERM
] = 1;
2124 gs
.sig
[eglsTERM
] = -1;
2129 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2130 signal_name
[last_signal_number_recvd
],
2131 gs
.sig
[eglsTERM
]==-1 ? "NS " : "");
2135 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2136 signal_name
[last_signal_number_recvd
],
2137 gs
.sig
[eglsTERM
]==-1 ? "NS " : "");
2139 handledSignal
=last_signal_number_recvd
;
2141 else if (MASTER(cr
) && (bNS
|| ir
->nstlist
<= 0) &&
2142 (max_hours
> 0 && run_time
> max_hours
*60.0*60.0*0.99) &&
2143 gs
.sig
[eglsTERM
] == 0 && gs
.set
[eglsTERM
] == 0)
2145 /* Signal to terminate the run */
2146 gs
.sig
[eglsTERM
] = (ir
->nstlist
== 0 ? 1 : -1);
2149 fprintf(fplog
,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step
,sbuf
),max_hours
*0.99);
2151 fprintf(stderr
, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step
,sbuf
),max_hours
*0.99);
2154 if (bResetCountersHalfMaxH
&& MASTER(cr
) &&
2155 run_time
> max_hours
*60.0*60.0*0.495)
2157 gs
.sig
[eglsRESETCOUNTERS
] = 1;
2160 if (ir
->nstlist
== -1 && !bRerunMD
)
2162 /* When bGStatEveryStep=FALSE, global_stat is only called
2163 * when we check the atom displacements, not at NS steps.
2164 * This means that also the bonded interaction count check is not
2165 * performed immediately after NS. Therefore a few MD steps could
2166 * be performed with missing interactions.
2167 * But wrong energies are never written to file,
2168 * since energies are only written after global_stat
2171 if (step
>= nlh
.step_nscheck
)
2173 nlh
.nabnsb
= natoms_beyond_ns_buffer(ir
,fr
,&top
->cgs
,
2174 nlh
.scale_tot
,state
->x
);
2178 /* This is not necessarily true,
2179 * but step_nscheck is determined quite conservatively.
2185 /* In parallel we only have to check for checkpointing in steps
2186 * where we do global communication,
2187 * otherwise the other nodes don't know.
2189 if (MASTER(cr
) && ((bGStat
|| !PAR(cr
)) &&
2192 run_time
>= nchkpt
*cpt_period
*60.0)) &&
2193 gs
.set
[eglsCHKPT
] == 0)
2195 gs
.sig
[eglsCHKPT
] = 1;
2200 gmx_iterate_init(&iterate
,bIterations
);
2203 /* for iterations, we save these vectors, as we will be redoing the calculations */
2204 if (bIterations
&& iterate
.bIterate
)
2206 copy_coupling_state(state
,bufstate
,ekind
,ekind_save
,&(ir
->opts
));
2208 bFirstIterate
= TRUE
;
2209 while (bFirstIterate
|| (bIterations
&& iterate
.bIterate
))
2211 /* We now restore these vectors to redo the calculation with improved extended variables */
2214 copy_coupling_state(bufstate
,state
,ekind_save
,ekind
,&(ir
->opts
));
2217 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
2218 so scroll down for that logic */
2220 /* ######### START SECOND UPDATE STEP ################# */
2221 GMX_MPE_LOG(ev_update_start
);
2223 if (!bRerunMD
|| rerun_fr
.bV
|| bForceUpdate
)
2225 wallcycle_start(wcycle
,ewcUPDATE
);
2227 /* Box is changed in update() when we do pressure coupling,
2228 * but we should still use the old box for energy corrections and when
2229 * writing it to the energy file, so it matches the trajectory files for
2230 * the same timestep above. Make a copy in a separate array.
2232 copy_mat(state
->box
,lastbox
);
2233 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
2236 if (bIterations
&& iterate
.bIterate
)
2244 /* we use a new value of scalevir to converge the iterations faster */
2245 scalevir
= tracevir
/trace(shake_vir
);
2247 msmul(shake_vir
,scalevir
,shake_vir
);
2248 m_add(force_vir
,shake_vir
,total_vir
);
2249 clear_mat(shake_vir
);
2251 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[3]);
2253 /* We can only do Berendsen coupling after we have summed
2254 * the kinetic energy or virial. Since the happens
2255 * in global_state after update, we should only do it at
2256 * step % nstlist = 1 with bGStatEveryStep=FALSE.
2259 update_extended(fplog
,step
,ir
,mdatoms
,state
,ekind
,pcoupl_mu
,M
,wcycle
,
2260 upd
,bInitStep
,FALSE
,&MassQ
);
2262 /* velocity (for VV) */
2263 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2264 ekind
,M
,wcycle
,upd
,FALSE
,etrtVELOCITY
,cr
,nrnb
,constr
,&top
->idef
);
2266 /* Above, initialize just copies ekinh into ekin,
2267 * it doesn't copy position (for VV),
2268 * and entire integrator for MD.
2273 copy_rvecn(state
->x
,cbuf
,0,state
->natoms
);
2276 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2277 ekind
,M
,wcycle
,upd
,bInitStep
,etrtPOSITION
,cr
,nrnb
,constr
,&top
->idef
);
2278 wallcycle_stop(wcycle
,ewcUPDATE
);
2280 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
2281 &top
->idef
,shake_vir
,force_vir
,
2282 cr
,nrnb
,wcycle
,upd
,constr
,
2283 bInitStep
,FALSE
,bCalcEnerPres
,state
->veta
);
2287 /* erase F_EKIN and F_TEMP here? */
2288 /* just compute the kinetic energy at the half step to perform a trotter step */
2289 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2290 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
2291 constr
,NULL
,FALSE
,lastbox
,
2292 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2293 cglo_flags
| CGLO_TEMPERATURE
| CGLO_CONSTRAINT
2295 wallcycle_start(wcycle
,ewcUPDATE
);
2296 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[4]);
2297 /* now we know the scaling, we can compute the positions again again */
2298 copy_rvecn(cbuf
,state
->x
,0,state
->natoms
);
2300 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2301 ekind
,M
,wcycle
,upd
,bInitStep
,etrtPOSITION
,cr
,nrnb
,constr
,&top
->idef
);
2302 wallcycle_stop(wcycle
,ewcUPDATE
);
2304 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
2305 /* are the small terms in the shake_vir here due
2306 * to numerical errors, or are they important
2307 * physically? I'm thinking they are just errors, but not completely sure.
2308 * For now, will call without actually constraining, constr=NULL*/
2309 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
2310 &top
->idef
,tmp_vir
,force_vir
,
2311 cr
,nrnb
,wcycle
,upd
,NULL
,
2312 bInitStep
,FALSE
,bCalcEnerPres
,
2315 if (!bOK
&& !bFFscan
)
2317 gmx_fatal(FARGS
,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2320 if (fr
->bSepDVDL
&& fplog
&& do_log
)
2322 fprintf(fplog
,sepdvdlformat
,"Constraint",0.0,dvdl
);
2324 enerd
->term
[F_DHDL_CON
] += dvdl
;
2328 /* Need to unshift here */
2329 unshift_self(graph
,state
->box
,state
->x
);
2332 GMX_BARRIER(cr
->mpi_comm_mygroup
);
2333 GMX_MPE_LOG(ev_update_finish
);
2337 wallcycle_start(wcycle
,ewcVSITECONSTR
);
2340 shift_self(graph
,state
->box
,state
->x
);
2342 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,state
->v
,
2343 top
->idef
.iparams
,top
->idef
.il
,
2344 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
2348 unshift_self(graph
,state
->box
,state
->x
);
2350 wallcycle_stop(wcycle
,ewcVSITECONSTR
);
2353 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
2354 if (ir
->nstlist
== -1 && bFirstIterate
)
2356 gs
.sig
[eglsNABNSB
] = nlh
.nabnsb
;
2358 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2359 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
2361 bFirstIterate
? &gs
: NULL
,(step
% gs
.nstms
== 0),
2363 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2365 | (!EI_VV(ir
->eI
) ? CGLO_ENERGY
: 0)
2366 | (!EI_VV(ir
->eI
) ? CGLO_TEMPERATURE
: 0)
2367 | (!EI_VV(ir
->eI
) || bRerunMD
? CGLO_PRESSURE
: 0)
2368 | (bIterations
&& iterate
.bIterate
? CGLO_ITERATE
: 0)
2369 | (bFirstIterate
? CGLO_FIRSTITERATE
: 0)
2372 if (ir
->nstlist
== -1 && bFirstIterate
)
2374 nlh
.nabnsb
= gs
.set
[eglsNABNSB
];
2375 gs
.set
[eglsNABNSB
] = 0;
2377 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
2378 /* ############# END CALC EKIN AND PRESSURE ################# */
2380 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
2381 the virial that should probably be addressed eventually. state->veta has better properies,
2382 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
2383 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
2386 done_iterating(cr
,fplog
,step
,&iterate
,bFirstIterate
,
2387 trace(shake_vir
),&tracevir
))
2391 bFirstIterate
= FALSE
;
2394 update_box(fplog
,step
,ir
,mdatoms
,state
,graph
,f
,
2395 ir
->nstlist
==-1 ? &nlh
.scale_tot
: NULL
,pcoupl_mu
,nrnb
,wcycle
,upd
,bInitStep
,FALSE
);
2397 /* ################# END UPDATE STEP 2 ################# */
2398 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
2400 /* The coordinates (x) were unshifted in update */
2401 if (bFFscan
&& (shellfc
==NULL
|| bConverged
))
2403 if (print_forcefield(fplog
,enerd
->term
,mdatoms
->homenr
,
2405 &(top_global
->mols
),mdatoms
->massT
,pres
))
2407 if (gmx_parallel_env_initialized())
2411 fprintf(stderr
,"\n");
2417 /* We will not sum ekinh_old,
2418 * so signal that we still have to do it.
2420 bSumEkinhOld
= TRUE
;
2425 /* Only do GCT when the relaxation of shells (minimization) has converged,
2426 * otherwise we might be coupling to bogus energies.
2427 * In parallel we must always do this, because the other sims might
2431 /* Since this is called with the new coordinates state->x, I assume
2432 * we want the new box state->box too. / EL 20040121
2434 do_coupling(fplog
,oenv
,nfile
,fnm
,tcr
,t
,step
,enerd
->term
,fr
,
2436 mdatoms
,&(top
->idef
),mu_aver
,
2437 top_global
->mols
.nr
,cr
,
2438 state
->box
,total_vir
,pres
,
2439 mu_tot
,state
->x
,f
,bConverged
);
2443 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
2445 sum_dhdl(enerd
,state
->lambda
,ir
);
2446 /* use the directly determined last velocity, not actually the averaged half steps */
2447 if (bTrotter
&& ir
->eI
==eiVV
)
2449 enerd
->term
[F_EKIN
] = last_ekin
;
2451 enerd
->term
[F_ETOT
] = enerd
->term
[F_EPOT
] + enerd
->term
[F_EKIN
];
2460 if (IR_NVT_TROTTER(ir
)) {
2461 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] + last_conserved
;
2463 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] +
2464 NPT_energy(ir
,state
,&MassQ
);
2468 enerd
->term
[F_ECONSERVED
] =
2469 enerd
->term
[F_ETOT
] + vrescale_energy(&(ir
->opts
),
2470 state
->therm_integral
);
2476 /* Check for excessively large energies */
2480 real etot_max
= 1e200
;
2482 real etot_max
= 1e30
;
2484 if (fabs(enerd
->term
[F_ETOT
]) > etot_max
)
2486 fprintf(stderr
,"Energy too large (%g), giving up\n",
2487 enerd
->term
[F_ETOT
]);
2490 /* ######### END PREPARING EDR OUTPUT ########### */
2492 /* Time for performance */
2493 if (((step
% stepout
) == 0) || bLastStep
)
2495 runtime_upd_proc(runtime
);
2503 if (!(bStartingFromCpt
&& (EI_VV(ir
->eI
))))
2507 upd_mdebin(mdebin
,bDoDHDL
? outf
->fp_dhdl
: NULL
,TRUE
,
2508 t
,mdatoms
->tmass
,enerd
,state
,lastbox
,
2509 shake_vir
,force_vir
,total_vir
,pres
,
2510 ekind
,mu_tot
,constr
);
2514 upd_mdebin_step(mdebin
);
2517 do_dr
= do_per_step(step
,ir
->nstdisreout
);
2518 do_or
= do_per_step(step
,ir
->nstorireout
);
2520 print_ebin(outf
->fp_ene
,do_ene
,do_dr
,do_or
,do_log
?fplog
:NULL
,
2522 eprNORMAL
,bCompact
,mdebin
,fcd
,groups
,&(ir
->opts
));
2524 if (ir
->ePull
!= epullNO
)
2526 pull_print_output(ir
->pull
,step
,t
);
2529 if (do_per_step(step
,ir
->nstlog
))
2531 if(fflush(fplog
) != 0)
2533 gmx_fatal(FARGS
,"Cannot flush logfile - maybe you are out of quota?");
2539 /* Remaining runtime */
2540 if (MULTIMASTER(cr
) && do_verbose
)
2544 fprintf(stderr
,"\n");
2546 print_time(stderr
,runtime
,step
,ir
,cr
);
2549 /* Replica exchange */
2551 if ((repl_ex_nst
> 0) && (step
> 0) && !bLastStep
&&
2552 do_per_step(step
,repl_ex_nst
))
2554 bExchanged
= replica_exchange(fplog
,cr
,repl_ex
,
2555 state_global
,enerd
->term
,
2558 if (bExchanged
&& PAR(cr
))
2560 if (DOMAINDECOMP(cr
))
2562 dd_partition_system(fplog
,step
,cr
,TRUE
,1,
2563 state_global
,top_global
,ir
,
2564 state
,&f
,mdatoms
,top
,fr
,
2565 vsite
,shellfc
,constr
,
2570 bcast_state(cr
,state
,FALSE
);
2576 bStartingFromCpt
= FALSE
;
2578 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
2579 /* Complicated conditional when bGStatEveryStep=FALSE.
2580 * We can not just use bGStat, since then the simulation results
2581 * would depend on nstenergy and nstlog or step_nscheck.
2583 if (((state
->flags
& (1<<estPRES_PREV
)) ||
2584 (state
->flags
& (1<<estSVIR_PREV
)) ||
2585 (state
->flags
& (1<<estFVIR_PREV
))) &&
2587 (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0) ||
2588 (ir
->nstlist
< 0 && nlh
.nabnsb
> 0) ||
2589 (ir
->nstlist
== 0 && bGStat
)))
2591 /* Store the pressure in t_state for pressure coupling
2592 * at the next MD step.
2594 if (state
->flags
& (1<<estPRES_PREV
))
2596 copy_mat(pres
,state
->pres_prev
);
2600 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
2604 /* read next frame from input trajectory */
2605 bNotLastFrame
= read_next_frame(oenv
,status
,&rerun_fr
);
2608 if (!bRerunMD
|| !rerun_fr
.bStep
)
2610 /* increase the MD step number */
2615 cycles
= wallcycle_stop(wcycle
,ewcSTEP
);
2616 if (DOMAINDECOMP(cr
) && wcycle
)
2618 dd_cycles_add(cr
->dd
,cycles
,ddCyclStep
);
2621 if (step_rel
== wcycle_get_reset_counters(wcycle
) ||
2622 gs
.set
[eglsRESETCOUNTERS
] != 0)
2624 /* Reset all the counters related to performance over the run */
2625 reset_all_counters(fplog
,cr
,step
,&step_rel
,ir
,wcycle
,nrnb
,runtime
);
2626 wcycle_set_reset_counters(wcycle
,-1);
2627 bResetCountersHalfMaxH
= FALSE
;
2628 gs
.set
[eglsRESETCOUNTERS
] = 0;
2631 /* End of main MD loop */
2635 runtime_end(runtime
);
2642 if (!(cr
->duty
& DUTY_PME
))
2644 /* Tell the PME only node to finish */
2650 if (ir
->nstcalcenergy
> 0 && !bRerunMD
)
2652 print_ebin(outf
->fp_ene
,FALSE
,FALSE
,FALSE
,fplog
,step
,t
,
2653 eprAVER
,FALSE
,mdebin
,fcd
,groups
,&(ir
->opts
));
2661 if (ir
->nstlist
== -1 && nlh
.nns
> 0 && fplog
)
2663 fprintf(fplog
,"Average neighborlist lifetime: %.1f steps, std.dev.: %.1f steps\n",nlh
.s1
/nlh
.nns
,sqrt(nlh
.s2
/nlh
.nns
- sqr(nlh
.s1
/nlh
.nns
)));
2664 fprintf(fplog
,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh
.ab
/nlh
.nns
);
2667 if (shellfc
&& fplog
)
2669 fprintf(fplog
,"Fraction of iterations that converged: %.2f %%\n",
2670 (nconverged
*100.0)/step_rel
);
2671 fprintf(fplog
,"Average number of force evaluations per MD step: %.2f\n\n",
2675 if (repl_ex_nst
> 0 && MASTER(cr
))
2677 print_replica_exchange_statistics(fplog
,repl_ex
);
2680 runtime
->nsteps_done
= step_rel
;