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"
92 #include "gmx_membed.h"
102 #include "corewrap.h"
106 /* simulation conditions to transmit. Keep in mind that they are
107 transmitted to other nodes through an MPI_Reduce after
108 casting them to a real (so the signals can be sent together with other
109 data). This means that the only meaningful values are positive,
111 enum { eglsNABNSB
, eglsCHKPT
, eglsSTOPCOND
, eglsRESETCOUNTERS
, eglsNR
};
112 /* Is the signal in one simulation independent of other simulations? */
113 gmx_bool gs_simlocal
[eglsNR
] = { TRUE
, FALSE
, FALSE
, TRUE
};
116 int nstms
; /* The frequency for intersimulation communication */
117 int sig
[eglsNR
]; /* The signal set by one process in do_md */
118 int set
[eglsNR
]; /* The communicated signal, equal for all processes */
122 static int multisim_min(const gmx_multisim_t
*ms
,int nmin
,int n
)
125 gmx_bool bPos
,bEqual
;
130 gmx_sumi_sim(ms
->nsim
,buf
,ms
);
133 for(s
=0; s
<ms
->nsim
; s
++)
135 bPos
= bPos
&& (buf
[s
] > 0);
136 bEqual
= bEqual
&& (buf
[s
] == buf
[0]);
142 nmin
= min(nmin
,buf
[0]);
146 /* Find the least common multiple */
147 for(d
=2; d
<nmin
; d
++)
150 while (s
< ms
->nsim
&& d
% buf
[s
] == 0)
156 /* We found the LCM and it is less than nmin */
168 static int multisim_nstsimsync(const t_commrec
*cr
,
169 const t_inputrec
*ir
,int repl_ex_nst
)
176 nmin
= multisim_min(cr
->ms
,nmin
,ir
->nstlist
);
177 nmin
= multisim_min(cr
->ms
,nmin
,ir
->nstcalcenergy
);
178 nmin
= multisim_min(cr
->ms
,nmin
,repl_ex_nst
);
181 gmx_fatal(FARGS
,"Can not find an appropriate interval for inter-simulation communication, since nstlist, nstcalcenergy and -replex are all <= 0");
183 /* Avoid inter-simulation communication at every (second) step */
190 gmx_bcast(sizeof(int),&nmin
,cr
);
195 static void init_global_signals(globsig_t
*gs
,const t_commrec
*cr
,
196 const t_inputrec
*ir
,int repl_ex_nst
)
202 gs
->nstms
= multisim_nstsimsync(cr
,ir
,repl_ex_nst
);
205 fprintf(debug
,"Syncing simulations for checkpointing and termination every %d steps\n",gs
->nstms
);
213 for(i
=0; i
<eglsNR
; i
++)
220 static void copy_coupling_state(t_state
*statea
,t_state
*stateb
,
221 gmx_ekindata_t
*ekinda
,gmx_ekindata_t
*ekindb
, t_grpopts
* opts
)
224 /* MRS note -- might be able to get rid of some of the arguments. Look over it when it's all debugged */
228 /* Make sure we have enough space for x and v */
229 if (statea
->nalloc
> stateb
->nalloc
)
231 stateb
->nalloc
= statea
->nalloc
;
232 srenew(stateb
->x
,stateb
->nalloc
);
233 srenew(stateb
->v
,stateb
->nalloc
);
236 stateb
->natoms
= statea
->natoms
;
237 stateb
->ngtc
= statea
->ngtc
;
238 stateb
->nnhpres
= statea
->nnhpres
;
239 stateb
->veta
= statea
->veta
;
242 copy_mat(ekinda
->ekin
,ekindb
->ekin
);
243 for (i
=0; i
<stateb
->ngtc
; i
++)
245 ekindb
->tcstat
[i
].T
= ekinda
->tcstat
[i
].T
;
246 ekindb
->tcstat
[i
].Th
= ekinda
->tcstat
[i
].Th
;
247 copy_mat(ekinda
->tcstat
[i
].ekinh
,ekindb
->tcstat
[i
].ekinh
);
248 copy_mat(ekinda
->tcstat
[i
].ekinf
,ekindb
->tcstat
[i
].ekinf
);
249 ekindb
->tcstat
[i
].ekinscalef_nhc
= ekinda
->tcstat
[i
].ekinscalef_nhc
;
250 ekindb
->tcstat
[i
].ekinscaleh_nhc
= ekinda
->tcstat
[i
].ekinscaleh_nhc
;
251 ekindb
->tcstat
[i
].vscale_nhc
= ekinda
->tcstat
[i
].vscale_nhc
;
254 copy_rvecn(statea
->x
,stateb
->x
,0,stateb
->natoms
);
255 copy_rvecn(statea
->v
,stateb
->v
,0,stateb
->natoms
);
256 copy_mat(statea
->box
,stateb
->box
);
257 copy_mat(statea
->box_rel
,stateb
->box_rel
);
258 copy_mat(statea
->boxv
,stateb
->boxv
);
260 for (i
= 0; i
<stateb
->ngtc
; i
++)
262 nc
= i
*opts
->nhchainlength
;
263 for (j
=0; j
<opts
->nhchainlength
; j
++)
265 stateb
->nosehoover_xi
[nc
+j
] = statea
->nosehoover_xi
[nc
+j
];
266 stateb
->nosehoover_vxi
[nc
+j
] = statea
->nosehoover_vxi
[nc
+j
];
269 if (stateb
->nhpres_xi
!= NULL
)
271 for (i
= 0; i
<stateb
->nnhpres
; i
++)
273 nc
= i
*opts
->nhchainlength
;
274 for (j
=0; j
<opts
->nhchainlength
; j
++)
276 stateb
->nhpres_xi
[nc
+j
] = statea
->nhpres_xi
[nc
+j
];
277 stateb
->nhpres_vxi
[nc
+j
] = statea
->nhpres_vxi
[nc
+j
];
283 static real
compute_conserved_from_auxiliary(t_inputrec
*ir
, t_state
*state
, t_extmass
*MassQ
)
293 quantity
= NPT_energy(ir
,state
,MassQ
);
296 quantity
= vrescale_energy(&(ir
->opts
),state
->therm_integral
);
304 static void compute_globals(FILE *fplog
, gmx_global_stat_t gstat
, t_commrec
*cr
, t_inputrec
*ir
,
305 t_forcerec
*fr
, gmx_ekindata_t
*ekind
,
306 t_state
*state
, t_state
*state_global
, t_mdatoms
*mdatoms
,
307 t_nrnb
*nrnb
, t_vcm
*vcm
, gmx_wallcycle_t wcycle
,
308 gmx_enerdata_t
*enerd
,tensor force_vir
, tensor shake_vir
, tensor total_vir
,
309 tensor pres
, rvec mu_tot
, gmx_constr_t constr
,
310 globsig_t
*gs
,gmx_bool bInterSimGS
,
311 matrix box
, gmx_mtop_t
*top_global
, real
*pcurr
,
312 int natoms
, gmx_bool
*bSumEkinhOld
, int flags
)
316 tensor corr_vir
,corr_pres
,shakeall_vir
;
317 gmx_bool bEner
,bPres
,bTemp
, bVV
;
318 gmx_bool bRerunMD
, bStopCM
, bGStat
, bIterate
,
319 bFirstIterate
,bReadEkin
,bEkinAveVel
,bScaleEkin
, bConstrain
;
320 real ekin
,temp
,prescorr
,enercorr
,dvdlcorr
;
322 /* translate CGLO flags to gmx_booleans */
323 bRerunMD
= flags
& CGLO_RERUNMD
;
324 bStopCM
= flags
& CGLO_STOPCM
;
325 bGStat
= flags
& CGLO_GSTAT
;
327 bReadEkin
= (flags
& CGLO_READEKIN
);
328 bScaleEkin
= (flags
& CGLO_SCALEEKIN
);
329 bEner
= flags
& CGLO_ENERGY
;
330 bTemp
= flags
& CGLO_TEMPERATURE
;
331 bPres
= (flags
& CGLO_PRESSURE
);
332 bConstrain
= (flags
& CGLO_CONSTRAINT
);
333 bIterate
= (flags
& CGLO_ITERATE
);
334 bFirstIterate
= (flags
& CGLO_FIRSTITERATE
);
336 /* we calculate a full state kinetic energy either with full-step velocity verlet
337 or half step where we need the pressure */
339 bEkinAveVel
= (ir
->eI
==eiVV
|| (ir
->eI
==eiVVAK
&& bPres
) || bReadEkin
);
341 /* in initalization, it sums the shake virial in vv, and to
342 sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
344 /* ########## Kinetic energy ############## */
348 /* Non-equilibrium MD: this is parallellized, but only does communication
349 * when there really is NEMD.
352 if (PAR(cr
) && (ekind
->bNEMD
))
354 accumulate_u(cr
,&(ir
->opts
),ekind
);
359 restore_ekinstate_from_state(cr
,ekind
,&state_global
->ekinstate
);
364 calc_ke_part(state
,&(ir
->opts
),mdatoms
,ekind
,nrnb
,bEkinAveVel
,bIterate
);
369 /* Calculate center of mass velocity if necessary, also parallellized */
370 if (bStopCM
&& !bRerunMD
&& bEner
)
372 calc_vcm_grp(fplog
,mdatoms
->start
,mdatoms
->homenr
,mdatoms
,
373 state
->x
,state
->v
,vcm
);
377 if (bTemp
|| bPres
|| bEner
|| bConstrain
)
381 /* We will not sum ekinh_old,
382 * so signal that we still have to do it.
384 *bSumEkinhOld
= TRUE
;
391 for(i
=0; i
<eglsNR
; i
++)
393 gs_buf
[i
] = gs
->sig
[i
];
398 wallcycle_start(wcycle
,ewcMoveE
);
399 GMX_MPE_LOG(ev_global_stat_start
);
400 global_stat(fplog
,gstat
,cr
,enerd
,force_vir
,shake_vir
,mu_tot
,
402 gs
!= NULL
? eglsNR
: 0,gs_buf
,
404 *bSumEkinhOld
,flags
);
405 GMX_MPE_LOG(ev_global_stat_finish
);
406 wallcycle_stop(wcycle
,ewcMoveE
);
410 if (MULTISIM(cr
) && bInterSimGS
)
414 /* Communicate the signals between the simulations */
415 gmx_sum_sim(eglsNR
,gs_buf
,cr
->ms
);
417 /* Communicate the signals form the master to the others */
418 gmx_bcast(eglsNR
*sizeof(gs_buf
[0]),gs_buf
,cr
);
420 for(i
=0; i
<eglsNR
; i
++)
422 if (bInterSimGS
|| gs_simlocal
[i
])
424 /* Set the communicated signal only when it is non-zero,
425 * since signals might not be processed at each MD step.
427 gsi
= (gs_buf
[i
] >= 0 ?
428 (int)(gs_buf
[i
] + 0.5) :
429 (int)(gs_buf
[i
] - 0.5));
434 /* Turn off the local signal */
439 *bSumEkinhOld
= FALSE
;
443 if (!ekind
->bNEMD
&& debug
&& bTemp
&& (vcm
->nr
> 0))
446 mdatoms
->start
,mdatoms
->start
+mdatoms
->homenr
,
447 state
->v
,vcm
->group_p
[0],
448 mdatoms
->massT
,mdatoms
->tmass
,ekind
->ekin
);
452 /* Do center of mass motion removal */
453 if (bStopCM
&& !bRerunMD
) /* is this correct? Does it get called too often with this logic? */
455 check_cm_grp(fplog
,vcm
,ir
,1);
456 do_stopcm_grp(fplog
,mdatoms
->start
,mdatoms
->homenr
,mdatoms
->cVCM
,
457 state
->x
,state
->v
,vcm
);
458 inc_nrnb(nrnb
,eNR_STOPCM
,mdatoms
->homenr
);
464 /* Sum the kinetic energies of the groups & calc temp */
465 /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with IR_NPT_TROTTER */
466 /* three maincase: VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
467 Leap with AveVel is not supported; it's not clear that it will actually work.
468 bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
469 If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
470 bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
471 If FALSE, we go ahead and erase over it.
473 enerd
->term
[F_TEMP
] = sum_ekin(&(ir
->opts
),ekind
,&(enerd
->term
[F_DKDL
]),
474 bEkinAveVel
,bIterate
,bScaleEkin
);
476 enerd
->term
[F_EKIN
] = trace(ekind
->ekin
);
479 /* ########## Long range energy information ###### */
481 if (bEner
|| bPres
|| bConstrain
)
483 calc_dispcorr(fplog
,ir
,fr
,0,top_global
->natoms
,box
,state
->lambda
,
484 corr_pres
,corr_vir
,&prescorr
,&enercorr
,&dvdlcorr
);
487 if (bEner
&& bFirstIterate
)
489 enerd
->term
[F_DISPCORR
] = enercorr
;
490 enerd
->term
[F_EPOT
] += enercorr
;
491 enerd
->term
[F_DVDL
] += dvdlcorr
;
492 if (fr
->efep
!= efepNO
) {
493 enerd
->dvdl_lin
+= dvdlcorr
;
497 /* ########## Now pressure ############## */
498 if (bPres
|| bConstrain
)
501 m_add(force_vir
,shake_vir
,total_vir
);
503 /* Calculate pressure and apply LR correction if PPPM is used.
504 * Use the box from last timestep since we already called update().
507 enerd
->term
[F_PRES
] = calc_pres(fr
->ePBC
,ir
->nwall
,box
,ekind
->ekin
,total_vir
,pres
,
508 (fr
->eeltype
==eelPPPM
)?enerd
->term
[F_COUL_RECIP
]:0.0);
510 /* Calculate long range corrections to pressure and energy */
511 /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
512 and computes enerd->term[F_DISPCORR]. Also modifies the
513 total_vir and pres tesors */
515 m_add(total_vir
,corr_vir
,total_vir
);
516 m_add(pres
,corr_pres
,pres
);
517 enerd
->term
[F_PDISPCORR
] = prescorr
;
518 enerd
->term
[F_PRES
] += prescorr
;
519 *pcurr
= enerd
->term
[F_PRES
];
520 /* calculate temperature using virial */
521 enerd
->term
[F_VTEMP
] = calc_temp(trace(total_vir
),ir
->opts
.nrdf
[0]);
527 /* Definitions for convergence of iterated constraints */
529 /* iterate constraints up to 50 times */
530 #define MAXITERCONST 50
535 real f
,fprev
,x
,xprev
;
538 real allrelerr
[MAXITERCONST
+2];
539 int num_close
; /* number of "close" violations, caused by limited precision. */
543 #define CONVERGEITER 0.000000001
544 #define CLOSE_ENOUGH 0.000001000
546 #define CONVERGEITER 0.0001
547 #define CLOSE_ENOUGH 0.0050
550 /* we want to keep track of the close calls. If there are too many, there might be some other issues.
551 so we make sure that it's either less than some predetermined number, or if more than that number,
552 only some small fraction of the total. */
553 #define MAX_NUMBER_CLOSE 50
554 #define FRACTION_CLOSE 0.001
556 /* maximum length of cyclic traps to check, emerging from limited numerical precision */
559 static void gmx_iterate_init(gmx_iterate_t
*iterate
,gmx_bool bIterate
)
564 iterate
->bIterate
= bIterate
;
565 iterate
->num_close
= 0;
566 for (i
=0;i
<MAXITERCONST
+2;i
++)
568 iterate
->allrelerr
[i
] = 0;
572 static gmx_bool
done_iterating(const t_commrec
*cr
,FILE *fplog
, int nsteps
, gmx_iterate_t
*iterate
, gmx_bool bFirstIterate
, real fom
, real
*newf
)
574 /* monitor convergence, and use a secant search to propose new
577 The secant method computes x_{i+1} = x_{i} - f(x_{i}) * ---------------------
578 f(x_{i}) - f(x_{i-1})
580 The function we are trying to zero is fom-x, where fom is the
581 "figure of merit" which is the pressure (or the veta value) we
582 would get by putting in an old value of the pressure or veta into
583 the incrementor function for the step or half step. I have
584 verified that this gives the same answer as self consistent
585 iteration, usually in many fewer steps, especially for small tau_p.
587 We could possibly eliminate an iteration with proper use
588 of the value from the previous step, but that would take a bit
589 more bookkeeping, especially for veta, since tests indicate the
590 function of veta on the last step is not sufficiently close to
591 guarantee convergence this step. This is
592 good enough for now. On my tests, I could use tau_p down to
593 0.02, which is smaller that would ever be necessary in
594 practice. Generally, 3-5 iterations will be sufficient */
596 real relerr
,err
,xmin
;
604 iterate
->f
= fom
-iterate
->x
;
611 iterate
->f
= fom
-iterate
->x
; /* we want to zero this difference */
612 if ((iterate
->iter_i
> 1) && (iterate
->iter_i
< MAXITERCONST
))
614 if (iterate
->f
==iterate
->fprev
)
620 *newf
= iterate
->x
- (iterate
->x
-iterate
->xprev
)*(iterate
->f
)/(iterate
->f
-iterate
->fprev
);
625 /* just use self-consistent iteration the first step to initialize, or
626 if it's not converging (which happens occasionally -- need to investigate why) */
630 /* Consider a slight shortcut allowing us to exit one sooner -- we check the
631 difference between the closest of x and xprev to the new
632 value. To be 100% certain, we should check the difference between
633 the last result, and the previous result, or
635 relerr = (fabs((x-xprev)/fom));
637 but this is pretty much never necessary under typical conditions.
638 Checking numerically, it seems to lead to almost exactly the same
639 trajectories, but there are small differences out a few decimal
640 places in the pressure, and eventually in the v_eta, but it could
643 if (fabs(*newf-x) < fabs(*newf - xprev)) { xmin = x;} else { xmin = xprev;}
644 relerr = (fabs((*newf-xmin) / *newf));
647 err
= fabs((iterate
->f
-iterate
->fprev
));
648 relerr
= fabs(err
/fom
);
650 iterate
->allrelerr
[iterate
->iter_i
] = relerr
;
652 if (iterate
->iter_i
> 0)
656 fprintf(debug
,"Iterating NPT constraints: %6i %20.12f%14.6g%20.12f\n",
657 iterate
->iter_i
,fom
,relerr
,*newf
);
660 if ((relerr
< CONVERGEITER
) || (err
< CONVERGEITER
) || (fom
==0) || ((iterate
->x
== iterate
->xprev
) && iterate
->iter_i
> 1))
662 iterate
->bIterate
= FALSE
;
665 fprintf(debug
,"Iterating NPT constraints: CONVERGED\n");
669 if (iterate
->iter_i
> MAXITERCONST
)
671 if (relerr
< CLOSE_ENOUGH
)
674 for (i
=1;i
<CYCLEMAX
;i
++) {
675 if ((iterate
->allrelerr
[iterate
->iter_i
-(1+i
)] == iterate
->allrelerr
[iterate
->iter_i
-1]) &&
676 (iterate
->allrelerr
[iterate
->iter_i
-(1+i
)] == iterate
->allrelerr
[iterate
->iter_i
-(1+2*i
)])) {
680 fprintf(debug
,"Exiting from an NPT iterating cycle of length %d\n",i
);
687 /* step 1: trapped in a numerical attractor */
688 /* we are trapped in a numerical attractor, and can't converge any more, and are close to the final result.
689 Better to give up convergence here than have the simulation die.
691 iterate
->num_close
++;
696 /* Step #2: test if we are reasonably close for other reasons, then monitor the number. If not, die */
698 /* how many close calls have we had? If less than a few, we're OK */
699 if (iterate
->num_close
< MAX_NUMBER_CLOSE
)
701 sprintf(buf
,"Slight numerical convergence deviation with NPT at step %d, relative error only %10.5g, likely not a problem, continuing\n",nsteps
,relerr
);
702 md_print_warning(cr
,fplog
,buf
);
703 iterate
->num_close
++;
705 /* if more than a few, check the total fraction. If too high, die. */
706 } else if (iterate
->num_close
/(double)nsteps
> FRACTION_CLOSE
) {
707 gmx_fatal(FARGS
,"Could not converge NPT constraints, too many exceptions (%d%%\n",iterate
->num_close
/(double)nsteps
);
713 gmx_fatal(FARGS
,"Could not converge NPT constraints\n");
718 iterate
->xprev
= iterate
->x
;
720 iterate
->fprev
= iterate
->f
;
726 static void check_nst_param(FILE *fplog
,t_commrec
*cr
,
727 const char *desc_nst
,int nst
,
728 const char *desc_p
,int *p
)
732 if (*p
> 0 && *p
% nst
!= 0)
734 /* Round up to the next multiple of nst */
735 *p
= ((*p
)/nst
+ 1)*nst
;
736 sprintf(buf
,"NOTE: %s changes %s to %d\n",desc_nst
,desc_p
,*p
);
737 md_print_warning(cr
,fplog
,buf
);
741 static void reset_all_counters(FILE *fplog
,t_commrec
*cr
,
742 gmx_large_int_t step
,
743 gmx_large_int_t
*step_rel
,t_inputrec
*ir
,
744 gmx_wallcycle_t wcycle
,t_nrnb
*nrnb
,
745 gmx_runtime_t
*runtime
)
747 char buf
[STRLEN
],sbuf
[STEPSTRSIZE
];
749 /* Reset all the counters related to performance over the run */
750 sprintf(buf
,"Step %s: resetting all time and cycle counters\n",
751 gmx_step_str(step
,sbuf
));
752 md_print_warning(cr
,fplog
,buf
);
754 wallcycle_stop(wcycle
,ewcRUN
);
755 wallcycle_reset_all(wcycle
);
756 if (DOMAINDECOMP(cr
))
758 reset_dd_statistics_counters(cr
->dd
);
761 ir
->init_step
+= *step_rel
;
762 ir
->nsteps
-= *step_rel
;
764 wallcycle_start(wcycle
,ewcRUN
);
765 runtime_start(runtime
);
766 print_date_and_time(fplog
,cr
->nodeid
,"Restarted time",runtime
);
769 static void min_zero(int *n
,int i
)
771 if (i
> 0 && (*n
== 0 || i
< *n
))
777 static int lcd4(int i1
,int i2
,int i3
,int i4
)
788 gmx_incons("All 4 inputs for determininig nstglobalcomm are <= 0");
791 while (nst
> 1 && ((i1
> 0 && i1
% nst
!= 0) ||
792 (i2
> 0 && i2
% nst
!= 0) ||
793 (i3
> 0 && i3
% nst
!= 0) ||
794 (i4
> 0 && i4
% nst
!= 0)))
802 static int check_nstglobalcomm(FILE *fplog
,t_commrec
*cr
,
803 int nstglobalcomm
,t_inputrec
*ir
)
807 if (!EI_DYNAMICS(ir
->eI
))
812 if (nstglobalcomm
== -1)
814 if (!(ir
->nstcalcenergy
> 0 ||
820 if (ir
->nstenergy
> 0 && ir
->nstenergy
< nstglobalcomm
)
822 nstglobalcomm
= ir
->nstenergy
;
827 /* Ensure that we do timely global communication for
828 * (possibly) each of the four following options.
830 nstglobalcomm
= lcd4(ir
->nstcalcenergy
,
832 ir
->etc
!= etcNO
? ir
->nsttcouple
: 0,
833 ir
->epc
!= epcNO
? ir
->nstpcouple
: 0);
838 if (ir
->nstlist
> 0 &&
839 nstglobalcomm
> ir
->nstlist
&& nstglobalcomm
% ir
->nstlist
!= 0)
841 nstglobalcomm
= (nstglobalcomm
/ ir
->nstlist
)*ir
->nstlist
;
842 sprintf(buf
,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm
);
843 md_print_warning(cr
,fplog
,buf
);
845 if (ir
->nstcalcenergy
> 0)
847 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
848 "nstcalcenergy",&ir
->nstcalcenergy
);
850 if (ir
->etc
!= etcNO
&& ir
->nsttcouple
> 0)
852 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
853 "nsttcouple",&ir
->nsttcouple
);
855 if (ir
->epc
!= epcNO
&& ir
->nstpcouple
> 0)
857 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
858 "nstpcouple",&ir
->nstpcouple
);
861 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
862 "nstenergy",&ir
->nstenergy
);
864 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
865 "nstlog",&ir
->nstlog
);
868 if (ir
->comm_mode
!= ecmNO
&& ir
->nstcomm
< nstglobalcomm
)
870 sprintf(buf
,"WARNING: Changing nstcomm from %d to %d\n",
871 ir
->nstcomm
,nstglobalcomm
);
872 md_print_warning(cr
,fplog
,buf
);
873 ir
->nstcomm
= nstglobalcomm
;
876 return nstglobalcomm
;
879 void check_ir_old_tpx_versions(t_commrec
*cr
,FILE *fplog
,
880 t_inputrec
*ir
,gmx_mtop_t
*mtop
)
882 /* Check required for old tpx files */
883 if (IR_TWINRANGE(*ir
) && ir
->nstlist
> 1 &&
884 ir
->nstcalcenergy
% ir
->nstlist
!= 0)
886 md_print_warning(cr
,fplog
,"Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies");
888 if (gmx_mtop_ftype_count(mtop
,F_CONSTR
) +
889 gmx_mtop_ftype_count(mtop
,F_CONSTRNC
) > 0 &&
890 ir
->eConstrAlg
== econtSHAKE
)
892 md_print_warning(cr
,fplog
,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect");
893 if (ir
->epc
!= epcNO
)
895 gmx_fatal(FARGS
,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
898 check_nst_param(fplog
,cr
,"nstlist",ir
->nstlist
,
899 "nstcalcenergy",&ir
->nstcalcenergy
);
900 if (ir
->epc
!= epcNO
)
902 check_nst_param(fplog
,cr
,"nstlist",ir
->nstlist
,
903 "nstpcouple",&ir
->nstpcouple
);
905 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
906 "nstenergy",&ir
->nstenergy
);
907 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
908 "nstlog",&ir
->nstlog
);
909 if (ir
->efep
!= efepNO
)
911 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
912 "nstdhdl",&ir
->nstdhdl
);
918 gmx_bool bGStatEveryStep
;
919 gmx_large_int_t step_ns
;
920 gmx_large_int_t step_nscheck
;
931 static void reset_nlistheuristics(gmx_nlheur_t
*nlh
,gmx_large_int_t step
)
935 nlh
->step_nscheck
= step
;
938 static void init_nlistheuristics(gmx_nlheur_t
*nlh
,
939 gmx_bool bGStatEveryStep
,gmx_large_int_t step
)
941 nlh
->bGStatEveryStep
= bGStatEveryStep
;
948 reset_nlistheuristics(nlh
,step
);
951 static void update_nliststatistics(gmx_nlheur_t
*nlh
,gmx_large_int_t step
)
953 gmx_large_int_t nl_lt
;
954 char sbuf
[STEPSTRSIZE
],sbuf2
[STEPSTRSIZE
];
956 /* Determine the neighbor list life time */
957 nl_lt
= step
- nlh
->step_ns
;
960 fprintf(debug
,"%d atoms beyond ns buffer, updating neighbor list after %s steps\n",nlh
->nabnsb
,gmx_step_str(nl_lt
,sbuf
));
964 nlh
->s2
+= nl_lt
*nl_lt
;
965 nlh
->ab
+= nlh
->nabnsb
;
966 if (nlh
->lt_runav
== 0)
968 nlh
->lt_runav
= nl_lt
;
969 /* Initialize the fluctuation average
970 * such that at startup we check after 0 steps.
972 nlh
->lt_runav2
= sqr(nl_lt
/2.0);
974 /* Running average with 0.9 gives an exp. history of 9.5 */
975 nlh
->lt_runav2
= 0.9*nlh
->lt_runav2
+ 0.1*sqr(nlh
->lt_runav
- nl_lt
);
976 nlh
->lt_runav
= 0.9*nlh
->lt_runav
+ 0.1*nl_lt
;
977 if (nlh
->bGStatEveryStep
)
979 /* Always check the nlist validity */
980 nlh
->step_nscheck
= step
;
984 /* We check after: <life time> - 2*sigma
985 * The factor 2 is quite conservative,
986 * but we assume that with nstlist=-1 the user
987 * prefers exact integration over performance.
989 nlh
->step_nscheck
= step
990 + (int)(nlh
->lt_runav
- 2.0*sqrt(nlh
->lt_runav2
)) - 1;
994 fprintf(debug
,"nlist life time %s run av. %4.1f sig %3.1f check %s check with -gcom %d\n",
995 gmx_step_str(nl_lt
,sbuf
),nlh
->lt_runav
,sqrt(nlh
->lt_runav2
),
996 gmx_step_str(nlh
->step_nscheck
-step
+1,sbuf2
),
997 (int)(nlh
->lt_runav
- 2.0*sqrt(nlh
->lt_runav2
)));
1001 static void set_nlistheuristics(gmx_nlheur_t
*nlh
,gmx_bool bReset
,gmx_large_int_t step
)
1007 reset_nlistheuristics(nlh
,step
);
1011 update_nliststatistics(nlh
,step
);
1014 nlh
->step_ns
= step
;
1015 /* Initialize the cumulative coordinate scaling matrix */
1016 clear_mat(nlh
->scale_tot
);
1017 for(d
=0; d
<DIM
; d
++)
1019 nlh
->scale_tot
[d
][d
] = 1.0;
1023 static void rerun_parallel_comm(t_commrec
*cr
,t_trxframe
*fr
,
1024 gmx_bool
*bNotLastFrame
)
1029 bAlloc
= (fr
->natoms
== 0);
1031 if (MASTER(cr
) && !*bNotLastFrame
)
1037 gmx_bcast(sizeof(*fr
),fr
,cr
);
1041 *bNotLastFrame
= (fr
->natoms
>= 0);
1043 if (*bNotLastFrame
&& PARTDECOMP(cr
))
1045 /* x and v are the only variable size quantities stored in trr
1046 * that are required for rerun (f is not needed).
1050 snew(fr
->x
,fr
->natoms
);
1051 snew(fr
->v
,fr
->natoms
);
1055 gmx_bcast(fr
->natoms
*sizeof(fr
->x
[0]),fr
->x
[0],cr
);
1059 gmx_bcast(fr
->natoms
*sizeof(fr
->v
[0]),fr
->v
[0],cr
);
1064 double do_md(FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
1065 const output_env_t oenv
, gmx_bool bVerbose
,gmx_bool bCompact
,
1067 gmx_vsite_t
*vsite
,gmx_constr_t constr
,
1068 int stepout
,t_inputrec
*ir
,
1069 gmx_mtop_t
*top_global
,
1071 t_state
*state_global
,
1073 t_nrnb
*nrnb
,gmx_wallcycle_t wcycle
,
1074 gmx_edsam_t ed
,t_forcerec
*fr
,
1075 int repl_ex_nst
,int repl_ex_seed
,gmx_membed_t
*membed
,
1076 real cpt_period
,real max_hours
,
1077 const char *deviceOptions
,
1078 unsigned long Flags
,
1079 gmx_runtime_t
*runtime
)
1082 gmx_large_int_t step
,step_rel
;
1085 gmx_bool bGStatEveryStep
,bGStat
,bNstEner
,bCalcEnerPres
;
1086 gmx_bool bNS
,bNStList
,bSimAnn
,bStopCM
,bRerunMD
,bNotLastFrame
=FALSE
,
1087 bFirstStep
,bStateFromTPX
,bInitStep
,bLastStep
,
1088 bBornRadii
,bStartingFromCpt
;
1089 gmx_bool bDoDHDL
=FALSE
;
1090 gmx_bool do_ene
,do_log
,do_verbose
,bRerunWarnNoV
=TRUE
,
1091 bForceUpdate
=FALSE
,bCPT
;
1093 gmx_bool bMasterState
;
1094 int force_flags
,cglo_flags
;
1095 tensor force_vir
,shake_vir
,total_vir
,tmp_vir
,pres
;
1097 t_trxstatus
*status
;
1100 t_state
*bufstate
=NULL
;
1101 matrix
*scale_tot
,pcoupl_mu
,M
,ebox
;
1103 t_trxframe rerun_fr
;
1104 gmx_repl_ex_t repl_ex
=NULL
;
1107 gmx_localtop_t
*top
;
1108 t_mdebin
*mdebin
=NULL
;
1109 t_state
*state
=NULL
;
1110 rvec
*f_global
=NULL
;
1113 gmx_enerdata_t
*enerd
;
1115 gmx_global_stat_t gstat
;
1116 gmx_update_t upd
=NULL
;
1117 t_graph
*graph
=NULL
;
1121 gmx_groups_t
*groups
;
1122 gmx_ekindata_t
*ekind
, *ekind_save
;
1123 gmx_shellfc_t shellfc
;
1124 int count
,nconverged
=0;
1127 gmx_bool bIonize
=FALSE
;
1128 gmx_bool bTCR
=FALSE
,bConverged
=TRUE
,bOK
,bSumEkinhOld
,bExchanged
;
1130 gmx_bool bResetCountersHalfMaxH
=FALSE
;
1131 gmx_bool bVV
,bIterations
,bFirstIterate
,bTemp
,bPres
,bTrotter
;
1132 real temp0
,mu_aver
=0,dvdl
;
1134 atom_id
*grpindex
=NULL
;
1136 t_coupl_rec
*tcr
=NULL
;
1137 rvec
*xcopy
=NULL
,*vcopy
=NULL
,*cbuf
=NULL
;
1138 matrix boxcopy
={{0}},lastbox
;
1140 real fom
,oldfom
,veta_save
,pcurr
,scalevir
,tracevir
;
1143 real saved_conserved_quantity
= 0;
1148 char sbuf
[STEPSTRSIZE
],sbuf2
[STEPSTRSIZE
];
1149 int handled_stop_condition
=gmx_stop_cond_none
; /* compare to get_stop_condition*/
1150 gmx_iterate_t iterate
;
1152 /* Temporary addition for FAHCORE checkpointing */
1156 /* Check for special mdrun options */
1157 bRerunMD
= (Flags
& MD_RERUN
);
1158 bIonize
= (Flags
& MD_IONIZE
);
1159 bFFscan
= (Flags
& MD_FFSCAN
);
1160 bAppend
= (Flags
& MD_APPENDFILES
);
1161 if (Flags
& MD_RESETCOUNTERSHALFWAY
)
1165 /* Signal to reset the counters half the simulation steps. */
1166 wcycle_set_reset_counters(wcycle
,ir
->nsteps
/2);
1168 /* Signal to reset the counters halfway the simulation time. */
1169 bResetCountersHalfMaxH
= (max_hours
> 0);
1172 /* md-vv uses averaged full step velocities for T-control
1173 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
1174 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
1175 bVV
= EI_VV(ir
->eI
);
1176 if (bVV
) /* to store the initial velocities while computing virial */
1178 snew(cbuf
,top_global
->natoms
);
1180 /* all the iteratative cases - only if there are constraints */
1181 bIterations
= ((IR_NPT_TROTTER(ir
)) && (constr
) && (!bRerunMD
));
1182 bTrotter
= (bVV
&& (IR_NPT_TROTTER(ir
) || (IR_NVT_TROTTER(ir
))));
1186 /* Since we don't know if the frames read are related in any way,
1187 * rebuild the neighborlist at every step.
1190 ir
->nstcalcenergy
= 1;
1194 check_ir_old_tpx_versions(cr
,fplog
,ir
,top_global
);
1196 nstglobalcomm
= check_nstglobalcomm(fplog
,cr
,nstglobalcomm
,ir
);
1197 bGStatEveryStep
= (nstglobalcomm
== 1);
1199 if (!bGStatEveryStep
&& ir
->nstlist
== -1 && fplog
!= NULL
)
1202 "To reduce the energy communication with nstlist = -1\n"
1203 "the neighbor list validity should not be checked at every step,\n"
1204 "this means that exact integration is not guaranteed.\n"
1205 "The neighbor list validity is checked after:\n"
1206 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
1207 "In most cases this will result in exact integration.\n"
1208 "This reduces the energy communication by a factor of 2 to 3.\n"
1209 "If you want less energy communication, set nstlist > 3.\n\n");
1212 if (bRerunMD
|| bFFscan
)
1216 groups
= &top_global
->groups
;
1218 /* Initial values */
1219 init_md(fplog
,cr
,ir
,oenv
,&t
,&t0
,&state_global
->lambda
,&lam0
,
1220 nrnb
,top_global
,&upd
,
1221 nfile
,fnm
,&outf
,&mdebin
,
1222 force_vir
,shake_vir
,mu_tot
,&bSimAnn
,&vcm
,state_global
,Flags
);
1224 clear_mat(total_vir
);
1226 /* Energy terms and groups */
1228 init_enerdata(top_global
->groups
.grps
[egcENER
].nr
,ir
->n_flambda
,enerd
);
1229 if (DOMAINDECOMP(cr
))
1235 snew(f
,top_global
->natoms
);
1238 /* Kinetic energy data */
1240 init_ekindata(fplog
,top_global
,&(ir
->opts
),ekind
);
1241 /* needed for iteration of constraints */
1243 init_ekindata(fplog
,top_global
,&(ir
->opts
),ekind_save
);
1244 /* Copy the cos acceleration to the groups struct */
1245 ekind
->cosacc
.cos_accel
= ir
->cos_accel
;
1247 gstat
= global_stat_init(ir
);
1250 /* Check for polarizable models and flexible constraints */
1251 shellfc
= init_shell_flexcon(fplog
,
1252 top_global
,n_flexible_constraints(constr
),
1253 (ir
->bContinuation
||
1254 (DOMAINDECOMP(cr
) && !MASTER(cr
))) ?
1255 NULL
: state_global
->x
);
1260 tMPI_Thread_mutex_lock(&deform_init_box_mutex
);
1262 set_deform_reference_box(upd
,
1263 deform_init_init_step_tpx
,
1264 deform_init_box_tpx
);
1266 tMPI_Thread_mutex_unlock(&deform_init_box_mutex
);
1271 double io
= compute_io(ir
,top_global
->natoms
,groups
,mdebin
->ebin
->nener
,1);
1272 if ((io
> 2000) && MASTER(cr
))
1274 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
1278 if (DOMAINDECOMP(cr
)) {
1279 top
= dd_init_local_top(top_global
);
1282 dd_init_local_state(cr
->dd
,state_global
,state
);
1284 if (DDMASTER(cr
->dd
) && ir
->nstfout
) {
1285 snew(f_global
,state_global
->natoms
);
1289 /* Initialize the particle decomposition and split the topology */
1290 top
= split_system(fplog
,top_global
,ir
,cr
);
1292 pd_cg_range(cr
,&fr
->cg0
,&fr
->hcg
);
1293 pd_at_range(cr
,&a0
,&a1
);
1295 top
= gmx_mtop_generate_local_top(top_global
,ir
);
1298 a1
= top_global
->natoms
;
1301 state
= partdec_init_local_state(cr
,state_global
);
1304 atoms2md(top_global
,ir
,0,NULL
,a0
,a1
-a0
,mdatoms
);
1307 set_vsite_top(vsite
,top
,mdatoms
,cr
);
1310 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
) {
1311 graph
= mk_graph(fplog
,&(top
->idef
),0,top_global
->natoms
,FALSE
,FALSE
);
1315 make_local_shells(cr
,mdatoms
,shellfc
);
1318 if (ir
->pull
&& PAR(cr
)) {
1319 dd_make_local_pull_groups(NULL
,ir
->pull
,mdatoms
);
1323 if (DOMAINDECOMP(cr
))
1325 /* Distribute the charge groups over the nodes from the master node */
1326 dd_partition_system(fplog
,ir
->init_step
,cr
,TRUE
,1,
1327 state_global
,top_global
,ir
,
1328 state
,&f
,mdatoms
,top
,fr
,
1329 vsite
,shellfc
,constr
,
1333 update_mdatoms(mdatoms
,state
->lambda
);
1337 /* Update mdebin with energy history if appending to output files */
1338 if ( Flags
& MD_APPENDFILES
)
1340 restore_energyhistory_from_state(mdebin
,&state_global
->enerhist
);
1342 /* Set the initial energy history in state to zero by updating once */
1343 update_energyhistory(&state_global
->enerhist
,mdebin
);
1346 if ((state
->flags
& (1<<estLD_RNG
)) && (Flags
& MD_READ_RNG
)) {
1347 /* Set the random state if we read a checkpoint file */
1348 set_stochd_state(upd
,state
);
1351 /* Initialize constraints */
1353 if (!DOMAINDECOMP(cr
))
1354 set_constraints(constr
,top
,ir
,mdatoms
,cr
);
1357 /* Check whether we have to GCT stuff */
1358 bTCR
= ftp2bSet(efGCT
,nfile
,fnm
);
1361 fprintf(stderr
,"Will do General Coupling Theory!\n");
1363 gnx
= top_global
->mols
.nr
;
1365 for(i
=0; (i
<gnx
); i
++) {
1370 if (repl_ex_nst
> 0 && MASTER(cr
))
1371 repl_ex
= init_replica_exchange(fplog
,cr
->ms
,state_global
,ir
,
1372 repl_ex_nst
,repl_ex_seed
);
1374 if (!ir
->bContinuation
&& !bRerunMD
)
1376 if (mdatoms
->cFREEZE
&& (state
->flags
& (1<<estV
)))
1378 /* Set the velocities of frozen particles to zero */
1379 for(i
=mdatoms
->start
; i
<mdatoms
->start
+mdatoms
->homenr
; i
++)
1381 for(m
=0; m
<DIM
; m
++)
1383 if (ir
->opts
.nFreeze
[mdatoms
->cFREEZE
[i
]][m
])
1393 /* Constrain the initial coordinates and velocities */
1394 do_constrain_first(fplog
,constr
,ir
,mdatoms
,state
,f
,
1395 graph
,cr
,nrnb
,fr
,top
,shake_vir
);
1399 /* Construct the virtual sites for the initial configuration */
1400 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,NULL
,
1401 top
->idef
.iparams
,top
->idef
.il
,
1402 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
1408 /* I'm assuming we need global communication the first time! MRS */
1409 cglo_flags
= (CGLO_TEMPERATURE
| CGLO_GSTAT
1410 | (bVV
? CGLO_PRESSURE
:0)
1411 | (bVV
? CGLO_CONSTRAINT
:0)
1412 | (bRerunMD
? CGLO_RERUNMD
:0)
1413 | ((Flags
& MD_READ_EKIN
) ? CGLO_READEKIN
:0));
1415 bSumEkinhOld
= FALSE
;
1416 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
1417 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
1418 constr
,NULL
,FALSE
,state
->box
,
1419 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,cglo_flags
);
1420 if (ir
->eI
== eiVVAK
) {
1421 /* a second call to get the half step temperature initialized as well */
1422 /* we do the same call as above, but turn the pressure off -- internally to
1423 compute_globals, this is recognized as a velocity verlet half-step
1424 kinetic energy calculation. This minimized excess variables, but
1425 perhaps loses some logic?*/
1427 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
1428 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
1429 constr
,NULL
,FALSE
,state
->box
,
1430 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
1431 cglo_flags
&~ CGLO_PRESSURE
);
1434 /* Calculate the initial half step temperature, and save the ekinh_old */
1435 if (!(Flags
& MD_STARTFROMCPT
))
1437 for(i
=0; (i
<ir
->opts
.ngtc
); i
++)
1439 copy_mat(ekind
->tcstat
[i
].ekinh
,ekind
->tcstat
[i
].ekinh_old
);
1444 enerd
->term
[F_TEMP
] *= 2; /* result of averages being done over previous and current step,
1445 and there is no previous step */
1447 temp0
= enerd
->term
[F_TEMP
];
1449 /* if using an iterative algorithm, we need to create a working directory for the state. */
1452 bufstate
= init_bufstate(state
);
1456 snew(xcopy
,state
->natoms
);
1457 snew(vcopy
,state
->natoms
);
1458 copy_rvecn(state
->x
,xcopy
,0,state
->natoms
);
1459 copy_rvecn(state
->v
,vcopy
,0,state
->natoms
);
1460 copy_mat(state
->box
,boxcopy
);
1463 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
1464 temperature control */
1465 trotter_seq
= init_npt_vars(ir
,state
,&MassQ
,bTrotter
);
1469 if (constr
&& !ir
->bContinuation
&& ir
->eConstrAlg
== econtLINCS
)
1472 "RMS relative constraint deviation after constraining: %.2e\n",
1473 constr_rmsd(constr
,FALSE
));
1475 fprintf(fplog
,"Initial temperature: %g K\n",enerd
->term
[F_TEMP
]);
1478 fprintf(stderr
,"starting md rerun '%s', reading coordinates from"
1479 " input trajectory '%s'\n\n",
1480 *(top_global
->name
),opt2fn("-rerun",nfile
,fnm
));
1483 fprintf(stderr
,"Calculated time to finish depends on nsteps from "
1484 "run input file,\nwhich may not correspond to the time "
1485 "needed to process input trajectory.\n\n");
1491 fprintf(stderr
,"starting mdrun '%s'\n",
1492 *(top_global
->name
));
1493 if (ir
->nsteps
>= 0)
1495 sprintf(tbuf
,"%8.1f",(ir
->init_step
+ir
->nsteps
)*ir
->delta_t
);
1499 sprintf(tbuf
,"%s","infinite");
1501 if (ir
->init_step
> 0)
1503 fprintf(stderr
,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
1504 gmx_step_str(ir
->init_step
+ir
->nsteps
,sbuf
),tbuf
,
1505 gmx_step_str(ir
->init_step
,sbuf2
),
1506 ir
->init_step
*ir
->delta_t
);
1510 fprintf(stderr
,"%s steps, %s ps.\n",
1511 gmx_step_str(ir
->nsteps
,sbuf
),tbuf
);
1514 fprintf(fplog
,"\n");
1517 /* Set and write start time */
1518 runtime_start(runtime
);
1519 print_date_and_time(fplog
,cr
->nodeid
,"Started mdrun",runtime
);
1520 wallcycle_start(wcycle
,ewcRUN
);
1522 fprintf(fplog
,"\n");
1524 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
1526 chkpt_ret
=fcCheckPointParallel( cr
->nodeid
,
1528 if ( chkpt_ret
== 0 )
1529 gmx_fatal( 3,__FILE__
,__LINE__
, "Checkpoint error on step %d\n", 0 );
1533 /***********************************************************
1535 * Loop over MD steps
1537 ************************************************************/
1539 /* if rerunMD then read coordinates and velocities from input trajectory */
1542 if (getenv("GMX_FORCE_UPDATE"))
1544 bForceUpdate
= TRUE
;
1547 rerun_fr
.natoms
= 0;
1550 bNotLastFrame
= read_first_frame(oenv
,&status
,
1551 opt2fn("-rerun",nfile
,fnm
),
1552 &rerun_fr
,TRX_NEED_X
| TRX_READ_V
);
1553 if (rerun_fr
.natoms
!= top_global
->natoms
)
1556 "Number of atoms in trajectory (%d) does not match the "
1557 "run input file (%d)\n",
1558 rerun_fr
.natoms
,top_global
->natoms
);
1560 if (ir
->ePBC
!= epbcNONE
)
1564 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
);
1566 if (max_cutoff2(ir
->ePBC
,rerun_fr
.box
) < sqr(fr
->rlistlong
))
1568 gmx_fatal(FARGS
,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr
.step
,rerun_fr
.time
);
1575 rerun_parallel_comm(cr
,&rerun_fr
,&bNotLastFrame
);
1578 if (ir
->ePBC
!= epbcNONE
)
1580 /* Set the shift vectors.
1581 * Necessary here when have a static box different from the tpr box.
1583 calc_shifts(rerun_fr
.box
,fr
->shift_vec
);
1587 /* loop over MD steps or if rerunMD to end of input trajectory */
1589 /* Skip the first Nose-Hoover integration when we get the state from tpx */
1590 bStateFromTPX
= !opt2bSet("-cpi",nfile
,fnm
);
1591 bInitStep
= bFirstStep
&& (bStateFromTPX
|| bVV
);
1592 bStartingFromCpt
= (Flags
& MD_STARTFROMCPT
) && bInitStep
;
1594 bSumEkinhOld
= FALSE
;
1597 init_global_signals(&gs
,cr
,ir
,repl_ex_nst
);
1599 step
= ir
->init_step
;
1602 if (ir
->nstlist
== -1)
1604 init_nlistheuristics(&nlh
,bGStatEveryStep
,step
);
1607 bLastStep
= (bRerunMD
|| (ir
->nsteps
>= 0 && step_rel
> ir
->nsteps
));
1608 while (!bLastStep
|| (bRerunMD
&& bNotLastFrame
)) {
1610 wallcycle_start(wcycle
,ewcSTEP
);
1612 GMX_MPE_LOG(ev_timestep1
);
1615 if (rerun_fr
.bStep
) {
1616 step
= rerun_fr
.step
;
1617 step_rel
= step
- ir
->init_step
;
1619 if (rerun_fr
.bTime
) {
1629 bLastStep
= (step_rel
== ir
->nsteps
);
1630 t
= t0
+ step
*ir
->delta_t
;
1633 if (ir
->efep
!= efepNO
)
1635 if (bRerunMD
&& rerun_fr
.bLambda
&& (ir
->delta_lambda
!=0))
1637 state_global
->lambda
= rerun_fr
.lambda
;
1641 state_global
->lambda
= lam0
+ step
*ir
->delta_lambda
;
1643 state
->lambda
= state_global
->lambda
;
1644 bDoDHDL
= do_per_step(step
,ir
->nstdhdl
);
1649 update_annealing_target_temp(&(ir
->opts
),t
);
1654 if (!(DOMAINDECOMP(cr
) && !MASTER(cr
)))
1656 for(i
=0; i
<state_global
->natoms
; i
++)
1658 copy_rvec(rerun_fr
.x
[i
],state_global
->x
[i
]);
1662 for(i
=0; i
<state_global
->natoms
; i
++)
1664 copy_rvec(rerun_fr
.v
[i
],state_global
->v
[i
]);
1669 for(i
=0; i
<state_global
->natoms
; i
++)
1671 clear_rvec(state_global
->v
[i
]);
1675 fprintf(stderr
,"\nWARNING: Some frames do not contain velocities.\n"
1676 " Ekin, temperature and pressure are incorrect,\n"
1677 " the virial will be incorrect when constraints are present.\n"
1679 bRerunWarnNoV
= FALSE
;
1683 copy_mat(rerun_fr
.box
,state_global
->box
);
1684 copy_mat(state_global
->box
,state
->box
);
1686 if (vsite
&& (Flags
& MD_RERUN_VSITE
))
1688 if (DOMAINDECOMP(cr
))
1690 gmx_fatal(FARGS
,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
1694 /* Following is necessary because the graph may get out of sync
1695 * with the coordinates if we only have every N'th coordinate set
1697 mk_mshift(fplog
,graph
,fr
->ePBC
,state
->box
,state
->x
);
1698 shift_self(graph
,state
->box
,state
->x
);
1700 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,state
->v
,
1701 top
->idef
.iparams
,top
->idef
.il
,
1702 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
1705 unshift_self(graph
,state
->box
,state
->x
);
1710 /* Stop Center of Mass motion */
1711 bStopCM
= (ir
->comm_mode
!= ecmNO
&& do_per_step(step
,ir
->nstcomm
));
1713 /* Copy back starting coordinates in case we're doing a forcefield scan */
1716 for(ii
=0; (ii
<state
->natoms
); ii
++)
1718 copy_rvec(xcopy
[ii
],state
->x
[ii
]);
1719 copy_rvec(vcopy
[ii
],state
->v
[ii
]);
1721 copy_mat(boxcopy
,state
->box
);
1726 /* for rerun MD always do Neighbour Searching */
1727 bNS
= (bFirstStep
|| ir
->nstlist
!= 0);
1732 /* Determine whether or not to do Neighbour Searching and LR */
1733 bNStList
= (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0);
1735 bNS
= (bFirstStep
|| bExchanged
|| bNStList
||
1736 (ir
->nstlist
== -1 && nlh
.nabnsb
> 0));
1738 if (bNS
&& ir
->nstlist
== -1)
1740 set_nlistheuristics(&nlh
,bFirstStep
|| bExchanged
,step
);
1744 /* < 0 means stop at next step, > 0 means stop at next NS step */
1745 if ( (gs
.set
[eglsSTOPCOND
] < 0 ) ||
1746 ( (gs
.set
[eglsSTOPCOND
] > 0 ) && ( bNS
|| ir
->nstlist
==0)) )
1751 /* Determine whether or not to update the Born radii if doing GB */
1752 bBornRadii
=bFirstStep
;
1753 if (ir
->implicit_solvent
&& (step
% ir
->nstgbradii
==0))
1758 do_log
= do_per_step(step
,ir
->nstlog
) || bFirstStep
|| bLastStep
;
1759 do_verbose
= bVerbose
&&
1760 (step
% stepout
== 0 || bFirstStep
|| bLastStep
);
1762 if (bNS
&& !(bFirstStep
&& ir
->bContinuation
&& !bRerunMD
))
1766 bMasterState
= TRUE
;
1770 bMasterState
= FALSE
;
1771 /* Correct the new box if it is too skewed */
1772 if (DYNAMIC_BOX(*ir
))
1774 if (correct_box(fplog
,step
,state
->box
,graph
))
1776 bMasterState
= TRUE
;
1779 if (DOMAINDECOMP(cr
) && bMasterState
)
1781 dd_collect_state(cr
->dd
,state
,state_global
);
1785 if (DOMAINDECOMP(cr
))
1787 /* Repartition the domain decomposition */
1788 wallcycle_start(wcycle
,ewcDOMDEC
);
1789 dd_partition_system(fplog
,step
,cr
,
1790 bMasterState
,nstglobalcomm
,
1791 state_global
,top_global
,ir
,
1792 state
,&f
,mdatoms
,top
,fr
,
1793 vsite
,shellfc
,constr
,
1794 nrnb
,wcycle
,do_verbose
);
1795 wallcycle_stop(wcycle
,ewcDOMDEC
);
1796 /* If using an iterative integrator, reallocate space to match the decomposition */
1800 if (MASTER(cr
) && do_log
&& !bFFscan
)
1802 print_ebin_header(fplog
,step
,t
,state
->lambda
);
1805 if (ir
->efep
!= efepNO
)
1807 update_mdatoms(mdatoms
,state
->lambda
);
1810 if (bRerunMD
&& rerun_fr
.bV
)
1813 /* We need the kinetic energy at minus the half step for determining
1814 * the full step kinetic energy and possibly for T-coupling.*/
1815 /* This may not be quite working correctly yet . . . . */
1816 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
1817 wcycle
,enerd
,NULL
,NULL
,NULL
,NULL
,mu_tot
,
1818 constr
,NULL
,FALSE
,state
->box
,
1819 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
1820 CGLO_RERUNMD
| CGLO_GSTAT
| CGLO_TEMPERATURE
);
1822 clear_mat(force_vir
);
1824 /* Ionize the atoms if necessary */
1827 ionize(fplog
,oenv
,mdatoms
,top_global
,t
,ir
,state
->x
,state
->v
,
1828 mdatoms
->start
,mdatoms
->start
+mdatoms
->homenr
,state
->box
,cr
);
1831 /* Update force field in ffscan program */
1834 if (update_forcefield(fplog
,
1836 mdatoms
->nr
,state
->x
,state
->box
)) {
1837 if (gmx_parallel_env_initialized())
1845 GMX_MPE_LOG(ev_timestep2
);
1847 /* We write a checkpoint at this MD step when:
1848 * either at an NS step when we signalled through gs,
1849 * or at the last step (but not when we do not want confout),
1850 * but never at the first step or with rerun.
1852 bCPT
= (((gs
.set
[eglsCHKPT
] && bNS
) ||
1853 (bLastStep
&& (Flags
& MD_CONFOUT
))) &&
1854 step
> ir
->init_step
&& !bRerunMD
);
1857 gs
.set
[eglsCHKPT
] = 0;
1860 /* Determine the energy and pressure:
1861 * at nstcalcenergy steps and at energy output steps (set below).
1863 bNstEner
= do_per_step(step
,ir
->nstcalcenergy
);
1866 (ir
->epc
!= epcNO
&& do_per_step(step
,ir
->nstpcouple
)));
1868 /* Do we need global communication ? */
1869 bGStat
= (bCalcEnerPres
|| bStopCM
||
1870 do_per_step(step
,nstglobalcomm
) ||
1871 (ir
->nstlist
== -1 && !bRerunMD
&& step
>= nlh
.step_nscheck
));
1873 do_ene
= (do_per_step(step
,ir
->nstenergy
) || bLastStep
);
1875 if (do_ene
|| do_log
)
1877 bCalcEnerPres
= TRUE
;
1881 /* these CGLO_ options remain the same throughout the iteration */
1882 cglo_flags
= ((bRerunMD
? CGLO_RERUNMD
: 0) |
1883 (bStopCM
? CGLO_STOPCM
: 0) |
1884 (bGStat
? CGLO_GSTAT
: 0)
1887 force_flags
= (GMX_FORCE_STATECHANGED
|
1888 ((DYNAMIC_BOX(*ir
) || bRerunMD
) ? GMX_FORCE_DYNAMICBOX
: 0) |
1889 GMX_FORCE_ALLFORCES
|
1890 (bNStList
? GMX_FORCE_DOLR
: 0) |
1892 (bCalcEnerPres
? GMX_FORCE_VIRIAL
: 0) |
1893 (bDoDHDL
? GMX_FORCE_DHDL
: 0)
1898 /* Now is the time to relax the shells */
1899 count
=relax_shell_flexcon(fplog
,cr
,bVerbose
,bFFscan
? step
+1 : step
,
1901 bStopCM
,top
,top_global
,
1903 state
,f
,force_vir
,mdatoms
,
1904 nrnb
,wcycle
,graph
,groups
,
1905 shellfc
,fr
,bBornRadii
,t
,mu_tot
,
1906 state
->natoms
,&bConverged
,vsite
,
1917 /* The coordinates (x) are shifted (to get whole molecules)
1919 * This is parallellized as well, and does communication too.
1920 * Check comments in sim_util.c
1923 do_force(fplog
,cr
,ir
,step
,nrnb
,wcycle
,top
,top_global
,groups
,
1924 state
->box
,state
->x
,&state
->hist
,
1925 f
,force_vir
,mdatoms
,enerd
,fcd
,
1926 state
->lambda
,graph
,
1927 fr
,vsite
,mu_tot
,t
,outf
->fp_field
,ed
,bBornRadii
,
1928 (bNS
? GMX_FORCE_NS
: 0) | force_flags
);
1931 GMX_BARRIER(cr
->mpi_comm_mygroup
);
1935 mu_aver
= calc_mu_aver(cr
,state
->x
,mdatoms
->chargeA
,
1936 mu_tot
,&top_global
->mols
,mdatoms
,gnx
,grpindex
);
1939 if (bTCR
&& bFirstStep
)
1941 tcr
=init_coupling(fplog
,nfile
,fnm
,cr
,fr
,mdatoms
,&(top
->idef
));
1942 fprintf(fplog
,"Done init_coupling\n");
1946 if (bVV
&& !bStartingFromCpt
&& !bRerunMD
)
1947 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1949 if (ir
->eI
==eiVV
&& bInitStep
)
1951 /* if using velocity verlet with full time step Ekin,
1952 * take the first half step only to compute the
1953 * virial for the first step. From there,
1954 * revert back to the initial coordinates
1955 * so that the input is actually the initial step.
1957 copy_rvecn(state
->v
,cbuf
,0,state
->natoms
); /* should make this better for parallelizing? */
1959 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1960 trotter_update(ir
,step
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
,ettTSEQ1
);
1963 update_coords(fplog
,step
,ir
,mdatoms
,state
,
1964 f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
1965 ekind
,M
,wcycle
,upd
,bInitStep
,etrtVELOCITY1
,
1966 cr
,nrnb
,constr
,&top
->idef
);
1970 gmx_iterate_init(&iterate
,bIterations
&& !bInitStep
);
1972 /* for iterations, we save these vectors, as we will be self-consistently iterating
1975 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1977 /* save the state */
1978 if (bIterations
&& iterate
.bIterate
) {
1979 copy_coupling_state(state
,bufstate
,ekind
,ekind_save
,&(ir
->opts
));
1982 bFirstIterate
= TRUE
;
1983 while (bFirstIterate
|| (bIterations
&& iterate
.bIterate
))
1985 if (bIterations
&& iterate
.bIterate
)
1987 copy_coupling_state(bufstate
,state
,ekind_save
,ekind
,&(ir
->opts
));
1988 if (bFirstIterate
&& bTrotter
)
1990 /* The first time through, we need a decent first estimate
1991 of veta(t+dt) to compute the constraints. Do
1992 this by computing the box volume part of the
1993 trotter integration at this time. Nothing else
1994 should be changed by this routine here. If
1995 !(first time), we start with the previous value
1998 veta_save
= state
->veta
;
1999 trotter_update(ir
,step
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
,ettTSEQ0
);
2000 vetanew
= state
->veta
;
2001 state
->veta
= veta_save
;
2006 if ( !bRerunMD
|| rerun_fr
.bV
|| bForceUpdate
) { /* Why is rerun_fr.bV here? Unclear. */
2009 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
2010 &top
->idef
,shake_vir
,NULL
,
2011 cr
,nrnb
,wcycle
,upd
,constr
,
2012 bInitStep
,TRUE
,bCalcEnerPres
,vetanew
);
2014 if (!bOK
&& !bFFscan
)
2016 gmx_fatal(FARGS
,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2021 { /* Need to unshift here if a do_force has been
2022 called in the previous step */
2023 unshift_self(graph
,state
->box
,state
->x
);
2027 /* if VV, compute the pressure and constraints */
2028 /* For VV2, we strictly only need this if using pressure
2029 * control, but we really would like to have accurate pressures
2031 * Think about ways around this in the future?
2032 * For now, keep this choice in comments.
2034 /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
2035 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
2037 bTemp
= ((ir
->eI
==eiVV
&&(!bInitStep
)) || (ir
->eI
==eiVVAK
));
2038 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2039 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
2040 constr
,NULL
,FALSE
,state
->box
,
2041 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2044 | (bTemp
? CGLO_TEMPERATURE
:0)
2045 | (bPres
? CGLO_PRESSURE
: 0)
2046 | (bPres
? CGLO_CONSTRAINT
: 0)
2047 | ((bIterations
&& iterate
.bIterate
) ? CGLO_ITERATE
: 0)
2048 | (bFirstIterate
? CGLO_FIRSTITERATE
: 0)
2051 /* explanation of above:
2052 a) We compute Ekin at the full time step
2053 if 1) we are using the AveVel Ekin, and it's not the
2054 initial step, or 2) if we are using AveEkin, but need the full
2055 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
2056 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
2057 EkinAveVel because it's needed for the pressure */
2059 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
2064 trotter_update(ir
,step
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
,ettTSEQ2
);
2068 update_tcouple(fplog
,step
,ir
,state
,ekind
,wcycle
,upd
,&MassQ
,mdatoms
);
2073 done_iterating(cr
,fplog
,step
,&iterate
,bFirstIterate
,
2074 state
->veta
,&vetanew
))
2078 bFirstIterate
= FALSE
;
2081 if (bTrotter
&& !bInitStep
) {
2082 copy_mat(shake_vir
,state
->svir_prev
);
2083 copy_mat(force_vir
,state
->fvir_prev
);
2084 if (IR_NVT_TROTTER(ir
) && ir
->eI
==eiVV
) {
2085 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
2086 enerd
->term
[F_TEMP
] = sum_ekin(&(ir
->opts
),ekind
,NULL
,(ir
->eI
==eiVV
),FALSE
,FALSE
);
2087 enerd
->term
[F_EKIN
] = trace(ekind
->ekin
);
2090 /* if it's the initial step, we performed this first step just to get the constraint virial */
2091 if (bInitStep
&& ir
->eI
==eiVV
) {
2092 copy_rvecn(cbuf
,state
->v
,0,state
->natoms
);
2095 if (fr
->bSepDVDL
&& fplog
&& do_log
)
2097 fprintf(fplog
,sepdvdlformat
,"Constraint",0.0,dvdl
);
2099 enerd
->term
[F_DHDL_CON
] += dvdl
;
2101 GMX_MPE_LOG(ev_timestep1
);
2104 /* MRS -- now done iterating -- compute the conserved quantity */
2106 saved_conserved_quantity
= compute_conserved_from_auxiliary(ir
,state
,&MassQ
);
2109 last_ekin
= enerd
->term
[F_EKIN
]; /* does this get preserved through checkpointing? */
2111 if ((ir
->eDispCorr
!= edispcEnerPres
) && (ir
->eDispCorr
!= edispcAllEnerPres
))
2113 saved_conserved_quantity
-= enerd
->term
[F_DISPCORR
];
2117 /* ######## END FIRST UPDATE STEP ############## */
2118 /* ######## If doing VV, we now have v(dt) ###### */
2120 /* ################## START TRAJECTORY OUTPUT ################# */
2122 /* Now we have the energies and forces corresponding to the
2123 * coordinates at time t. We must output all of this before
2125 * for RerunMD t is read from input trajectory
2127 GMX_MPE_LOG(ev_output_start
);
2130 if (do_per_step(step
,ir
->nstxout
)) { mdof_flags
|= MDOF_X
; }
2131 if (do_per_step(step
,ir
->nstvout
)) { mdof_flags
|= MDOF_V
; }
2132 if (do_per_step(step
,ir
->nstfout
)) { mdof_flags
|= MDOF_F
; }
2133 if (do_per_step(step
,ir
->nstxtcout
)) { mdof_flags
|= MDOF_XTC
; }
2134 if (bCPT
) { mdof_flags
|= MDOF_CPT
; };
2136 #if defined(GMX_FAHCORE) || defined(GMX_WRITELASTSTEP)
2139 /* Enforce writing positions and velocities at end of run */
2140 mdof_flags
|= (MDOF_X
| MDOF_V
);
2145 fcReportProgress( ir
->nsteps
, step
);
2147 /* sync bCPT and fc record-keeping */
2148 if (bCPT
&& MASTER(cr
))
2149 fcRequestCheckPoint();
2152 if (mdof_flags
!= 0)
2154 wallcycle_start(wcycle
,ewcTRAJ
);
2157 if (state
->flags
& (1<<estLD_RNG
))
2159 get_stochd_state(upd
,state
);
2165 state_global
->ekinstate
.bUpToDate
= FALSE
;
2169 update_ekinstate(&state_global
->ekinstate
,ekind
);
2170 state_global
->ekinstate
.bUpToDate
= TRUE
;
2172 update_energyhistory(&state_global
->enerhist
,mdebin
);
2175 write_traj(fplog
,cr
,outf
,mdof_flags
,top_global
,
2176 step
,t
,state
,state_global
,f
,f_global
,&n_xtc
,&x_xtc
);
2183 if (bLastStep
&& step_rel
== ir
->nsteps
&&
2184 (Flags
& MD_CONFOUT
) && MASTER(cr
) &&
2185 !bRerunMD
&& !bFFscan
)
2187 /* x and v have been collected in write_traj,
2188 * because a checkpoint file will always be written
2191 fprintf(stderr
,"\nWriting final coordinates.\n");
2192 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
&&
2195 /* Make molecules whole only for confout writing */
2196 do_pbc_mtop(fplog
,ir
->ePBC
,state
->box
,top_global
,state_global
->x
);
2198 write_sto_conf_mtop(ftp2fn(efSTO
,nfile
,fnm
),
2199 *top_global
->name
,top_global
,
2200 state_global
->x
,state_global
->v
,
2201 ir
->ePBC
,state
->box
);
2204 wallcycle_stop(wcycle
,ewcTRAJ
);
2206 GMX_MPE_LOG(ev_output_finish
);
2208 /* kludge -- virial is lost with restart for NPT control. Must restart */
2209 if (bStartingFromCpt
&& bVV
)
2211 copy_mat(state
->svir_prev
,shake_vir
);
2212 copy_mat(state
->fvir_prev
,force_vir
);
2214 /* ################## END TRAJECTORY OUTPUT ################ */
2216 /* Determine the pressure:
2217 * always when we want exact averages in the energy file,
2218 * at ns steps when we have pressure coupling,
2219 * otherwise only at energy output steps (set below).
2222 bNstEner
= (bGStatEveryStep
|| do_per_step(step
,ir
->nstcalcenergy
));
2223 bCalcEnerPres
= bNstEner
;
2225 /* Do we need global communication ? */
2226 bGStat
= (bGStatEveryStep
|| bStopCM
|| bNS
||
2227 (ir
->nstlist
== -1 && !bRerunMD
&& step
>= nlh
.step_nscheck
));
2229 do_ene
= (do_per_step(step
,ir
->nstenergy
) || bLastStep
);
2231 if (do_ene
|| do_log
)
2233 bCalcEnerPres
= TRUE
;
2237 /* Determine the wallclock run time up till now */
2238 run_time
= gmx_gettime() - (double)runtime
->real
;
2240 /* Check whether everything is still allright */
2241 if (((int)gmx_get_stop_condition() > handled_stop_condition
)
2247 /* this is just make gs.sig compatible with the hack
2248 of sending signals around by MPI_Reduce with together with
2250 if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns
)
2251 gs
.sig
[eglsSTOPCOND
]=1;
2252 if ( gmx_get_stop_condition() == gmx_stop_cond_next
)
2253 gs
.sig
[eglsSTOPCOND
]=-1;
2254 /* < 0 means stop at next step, > 0 means stop at next NS step */
2258 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2259 gmx_get_signal_name(),
2260 gs
.sig
[eglsSTOPCOND
]==1 ? "NS " : "");
2264 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2265 gmx_get_signal_name(),
2266 gs
.sig
[eglsSTOPCOND
]==1 ? "NS " : "");
2268 handled_stop_condition
=(int)gmx_get_stop_condition();
2270 else if (MASTER(cr
) && (bNS
|| ir
->nstlist
<= 0) &&
2271 (max_hours
> 0 && run_time
> max_hours
*60.0*60.0*0.99) &&
2272 gs
.sig
[eglsSTOPCOND
] == 0 && gs
.set
[eglsSTOPCOND
] == 0)
2274 /* Signal to terminate the run */
2275 gs
.sig
[eglsSTOPCOND
] = 1;
2278 fprintf(fplog
,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step
,sbuf
),max_hours
*0.99);
2280 fprintf(stderr
, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step
,sbuf
),max_hours
*0.99);
2283 if (bResetCountersHalfMaxH
&& MASTER(cr
) &&
2284 run_time
> max_hours
*60.0*60.0*0.495)
2286 gs
.sig
[eglsRESETCOUNTERS
] = 1;
2289 if (ir
->nstlist
== -1 && !bRerunMD
)
2291 /* When bGStatEveryStep=FALSE, global_stat is only called
2292 * when we check the atom displacements, not at NS steps.
2293 * This means that also the bonded interaction count check is not
2294 * performed immediately after NS. Therefore a few MD steps could
2295 * be performed with missing interactions.
2296 * But wrong energies are never written to file,
2297 * since energies are only written after global_stat
2300 if (step
>= nlh
.step_nscheck
)
2302 nlh
.nabnsb
= natoms_beyond_ns_buffer(ir
,fr
,&top
->cgs
,
2303 nlh
.scale_tot
,state
->x
);
2307 /* This is not necessarily true,
2308 * but step_nscheck is determined quite conservatively.
2314 /* In parallel we only have to check for checkpointing in steps
2315 * where we do global communication,
2316 * otherwise the other nodes don't know.
2318 if (MASTER(cr
) && ((bGStat
|| !PAR(cr
)) &&
2321 run_time
>= nchkpt
*cpt_period
*60.0)) &&
2322 gs
.set
[eglsCHKPT
] == 0)
2324 gs
.sig
[eglsCHKPT
] = 1;
2329 gmx_iterate_init(&iterate
,bIterations
);
2332 /* for iterations, we save these vectors, as we will be redoing the calculations */
2333 if (bIterations
&& iterate
.bIterate
)
2335 copy_coupling_state(state
,bufstate
,ekind
,ekind_save
,&(ir
->opts
));
2337 bFirstIterate
= TRUE
;
2338 while (bFirstIterate
|| (bIterations
&& iterate
.bIterate
))
2340 /* We now restore these vectors to redo the calculation with improved extended variables */
2343 copy_coupling_state(bufstate
,state
,ekind_save
,ekind
,&(ir
->opts
));
2346 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
2347 so scroll down for that logic */
2349 /* ######### START SECOND UPDATE STEP ################# */
2350 GMX_MPE_LOG(ev_update_start
);
2351 /* Box is changed in update() when we do pressure coupling,
2352 * but we should still use the old box for energy corrections and when
2353 * writing it to the energy file, so it matches the trajectory files for
2354 * the same timestep above. Make a copy in a separate array.
2356 copy_mat(state
->box
,lastbox
);
2359 if (!(bRerunMD
&& !rerun_fr
.bV
&& !bForceUpdate
))
2361 wallcycle_start(wcycle
,ewcUPDATE
);
2363 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
2366 if (bIterations
&& iterate
.bIterate
)
2374 /* we use a new value of scalevir to converge the iterations faster */
2375 scalevir
= tracevir
/trace(shake_vir
);
2377 msmul(shake_vir
,scalevir
,shake_vir
);
2378 m_add(force_vir
,shake_vir
,total_vir
);
2379 clear_mat(shake_vir
);
2381 trotter_update(ir
,step
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
,ettTSEQ3
);
2382 /* We can only do Berendsen coupling after we have summed
2383 * the kinetic energy or virial. Since the happens
2384 * in global_state after update, we should only do it at
2385 * step % nstlist = 1 with bGStatEveryStep=FALSE.
2390 update_tcouple(fplog
,step
,ir
,state
,ekind
,wcycle
,upd
,&MassQ
,mdatoms
);
2391 update_pcouple(fplog
,step
,ir
,state
,pcoupl_mu
,M
,wcycle
,
2397 /* velocity half-step update */
2398 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,
2399 fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2400 ekind
,M
,wcycle
,upd
,FALSE
,etrtVELOCITY2
,
2401 cr
,nrnb
,constr
,&top
->idef
);
2404 /* Above, initialize just copies ekinh into ekin,
2405 * it doesn't copy position (for VV),
2406 * and entire integrator for MD.
2411 copy_rvecn(state
->x
,cbuf
,0,state
->natoms
);
2414 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2415 ekind
,M
,wcycle
,upd
,bInitStep
,etrtPOSITION
,cr
,nrnb
,constr
,&top
->idef
);
2416 wallcycle_stop(wcycle
,ewcUPDATE
);
2418 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
2419 &top
->idef
,shake_vir
,force_vir
,
2420 cr
,nrnb
,wcycle
,upd
,constr
,
2421 bInitStep
,FALSE
,bCalcEnerPres
,state
->veta
);
2425 /* erase F_EKIN and F_TEMP here? */
2426 /* just compute the kinetic energy at the half step to perform a trotter step */
2427 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2428 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
2429 constr
,NULL
,FALSE
,lastbox
,
2430 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2431 cglo_flags
| CGLO_TEMPERATURE
2433 wallcycle_start(wcycle
,ewcUPDATE
);
2434 trotter_update(ir
,step
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
,ettTSEQ4
);
2435 /* now we know the scaling, we can compute the positions again again */
2436 copy_rvecn(cbuf
,state
->x
,0,state
->natoms
);
2438 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2439 ekind
,M
,wcycle
,upd
,bInitStep
,etrtPOSITION
,cr
,nrnb
,constr
,&top
->idef
);
2440 wallcycle_stop(wcycle
,ewcUPDATE
);
2442 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
2443 /* are the small terms in the shake_vir here due
2444 * to numerical errors, or are they important
2445 * physically? I'm thinking they are just errors, but not completely sure.
2446 * For now, will call without actually constraining, constr=NULL*/
2447 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
2448 &top
->idef
,tmp_vir
,force_vir
,
2449 cr
,nrnb
,wcycle
,upd
,NULL
,
2450 bInitStep
,FALSE
,bCalcEnerPres
,
2453 if (!bOK
&& !bFFscan
)
2455 gmx_fatal(FARGS
,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2458 if (fr
->bSepDVDL
&& fplog
&& do_log
)
2460 fprintf(fplog
,sepdvdlformat
,"Constraint",0.0,dvdl
);
2462 enerd
->term
[F_DHDL_CON
] += dvdl
;
2466 /* Need to unshift here */
2467 unshift_self(graph
,state
->box
,state
->x
);
2470 GMX_BARRIER(cr
->mpi_comm_mygroup
);
2471 GMX_MPE_LOG(ev_update_finish
);
2475 wallcycle_start(wcycle
,ewcVSITECONSTR
);
2478 shift_self(graph
,state
->box
,state
->x
);
2480 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,state
->v
,
2481 top
->idef
.iparams
,top
->idef
.il
,
2482 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
2486 unshift_self(graph
,state
->box
,state
->x
);
2488 wallcycle_stop(wcycle
,ewcVSITECONSTR
);
2491 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
2492 if (ir
->nstlist
== -1 && bFirstIterate
)
2494 gs
.sig
[eglsNABNSB
] = nlh
.nabnsb
;
2496 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2497 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
2499 bFirstIterate
? &gs
: NULL
,(step
% gs
.nstms
== 0),
2501 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2503 | (!EI_VV(ir
->eI
) ? CGLO_ENERGY
: 0)
2504 | (!EI_VV(ir
->eI
) ? CGLO_TEMPERATURE
: 0)
2505 | (!EI_VV(ir
->eI
) || bRerunMD
? CGLO_PRESSURE
: 0)
2506 | (bIterations
&& iterate
.bIterate
? CGLO_ITERATE
: 0)
2507 | (bFirstIterate
? CGLO_FIRSTITERATE
: 0)
2510 if (ir
->nstlist
== -1 && bFirstIterate
)
2512 nlh
.nabnsb
= gs
.set
[eglsNABNSB
];
2513 gs
.set
[eglsNABNSB
] = 0;
2515 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
2516 /* ############# END CALC EKIN AND PRESSURE ################# */
2518 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
2519 the virial that should probably be addressed eventually. state->veta has better properies,
2520 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
2521 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
2524 done_iterating(cr
,fplog
,step
,&iterate
,bFirstIterate
,
2525 trace(shake_vir
),&tracevir
))
2529 bFirstIterate
= FALSE
;
2532 update_box(fplog
,step
,ir
,mdatoms
,state
,graph
,f
,
2533 ir
->nstlist
==-1 ? &nlh
.scale_tot
: NULL
,pcoupl_mu
,nrnb
,wcycle
,upd
,bInitStep
,FALSE
);
2535 /* ################# END UPDATE STEP 2 ################# */
2536 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
2538 /* The coordinates (x) were unshifted in update */
2539 if (bFFscan
&& (shellfc
==NULL
|| bConverged
))
2541 if (print_forcefield(fplog
,enerd
->term
,mdatoms
->homenr
,
2543 &(top_global
->mols
),mdatoms
->massT
,pres
))
2545 if (gmx_parallel_env_initialized())
2549 fprintf(stderr
,"\n");
2555 /* We will not sum ekinh_old,
2556 * so signal that we still have to do it.
2558 bSumEkinhOld
= TRUE
;
2563 /* Only do GCT when the relaxation of shells (minimization) has converged,
2564 * otherwise we might be coupling to bogus energies.
2565 * In parallel we must always do this, because the other sims might
2569 /* Since this is called with the new coordinates state->x, I assume
2570 * we want the new box state->box too. / EL 20040121
2572 do_coupling(fplog
,oenv
,nfile
,fnm
,tcr
,t
,step
,enerd
->term
,fr
,
2574 mdatoms
,&(top
->idef
),mu_aver
,
2575 top_global
->mols
.nr
,cr
,
2576 state
->box
,total_vir
,pres
,
2577 mu_tot
,state
->x
,f
,bConverged
);
2581 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
2583 /* sum up the foreign energy and dhdl terms */
2584 sum_dhdl(enerd
,state
->lambda
,ir
);
2586 /* use the directly determined last velocity, not actually the averaged half steps */
2587 if (bTrotter
&& ir
->eI
==eiVV
)
2589 enerd
->term
[F_EKIN
] = last_ekin
;
2591 enerd
->term
[F_ETOT
] = enerd
->term
[F_EPOT
] + enerd
->term
[F_EKIN
];
2595 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] + saved_conserved_quantity
;
2599 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] + compute_conserved_from_auxiliary(ir
,state
,&MassQ
);
2601 /* Check for excessively large energies */
2605 real etot_max
= 1e200
;
2607 real etot_max
= 1e30
;
2609 if (fabs(enerd
->term
[F_ETOT
]) > etot_max
)
2611 fprintf(stderr
,"Energy too large (%g), giving up\n",
2612 enerd
->term
[F_ETOT
]);
2615 /* ######### END PREPARING EDR OUTPUT ########### */
2617 /* Time for performance */
2618 if (((step
% stepout
) == 0) || bLastStep
)
2620 runtime_upd_proc(runtime
);
2626 gmx_bool do_dr
,do_or
;
2628 if (!(bStartingFromCpt
&& (EI_VV(ir
->eI
))))
2632 upd_mdebin(mdebin
,bDoDHDL
, TRUE
,
2633 t
,mdatoms
->tmass
,enerd
,state
,lastbox
,
2634 shake_vir
,force_vir
,total_vir
,pres
,
2635 ekind
,mu_tot
,constr
);
2639 upd_mdebin_step(mdebin
);
2642 do_dr
= do_per_step(step
,ir
->nstdisreout
);
2643 do_or
= do_per_step(step
,ir
->nstorireout
);
2645 print_ebin(outf
->fp_ene
,do_ene
,do_dr
,do_or
,do_log
?fplog
:NULL
,
2647 eprNORMAL
,bCompact
,mdebin
,fcd
,groups
,&(ir
->opts
));
2649 if (ir
->ePull
!= epullNO
)
2651 pull_print_output(ir
->pull
,step
,t
);
2654 if (do_per_step(step
,ir
->nstlog
))
2656 if(fflush(fplog
) != 0)
2658 gmx_fatal(FARGS
,"Cannot flush logfile - maybe you are out of quota?");
2664 /* Remaining runtime */
2665 if (MULTIMASTER(cr
) && (do_verbose
|| gmx_got_usr_signal() ))
2669 fprintf(stderr
,"\n");
2671 print_time(stderr
,runtime
,step
,ir
,cr
);
2674 /* Replica exchange */
2676 if ((repl_ex_nst
> 0) && (step
> 0) && !bLastStep
&&
2677 do_per_step(step
,repl_ex_nst
))
2679 bExchanged
= replica_exchange(fplog
,cr
,repl_ex
,
2680 state_global
,enerd
->term
,
2683 if (bExchanged
&& DOMAINDECOMP(cr
))
2685 dd_partition_system(fplog
,step
,cr
,TRUE
,1,
2686 state_global
,top_global
,ir
,
2687 state
,&f
,mdatoms
,top
,fr
,
2688 vsite
,shellfc
,constr
,
2695 bStartingFromCpt
= FALSE
;
2697 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
2698 /* Complicated conditional when bGStatEveryStep=FALSE.
2699 * We can not just use bGStat, since then the simulation results
2700 * would depend on nstenergy and nstlog or step_nscheck.
2702 if (((state
->flags
& (1<<estPRES_PREV
)) ||
2703 (state
->flags
& (1<<estSVIR_PREV
)) ||
2704 (state
->flags
& (1<<estFVIR_PREV
))) &&
2706 (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0) ||
2707 (ir
->nstlist
< 0 && nlh
.nabnsb
> 0) ||
2708 (ir
->nstlist
== 0 && bGStat
)))
2710 /* Store the pressure in t_state for pressure coupling
2711 * at the next MD step.
2713 if (state
->flags
& (1<<estPRES_PREV
))
2715 copy_mat(pres
,state
->pres_prev
);
2719 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
2721 if ( (membed
!=NULL
) && (!bLastStep
) )
2722 rescale_membed(step_rel
,membed
,state_global
->x
);
2728 /* read next frame from input trajectory */
2729 bNotLastFrame
= read_next_frame(oenv
,status
,&rerun_fr
);
2734 rerun_parallel_comm(cr
,&rerun_fr
,&bNotLastFrame
);
2738 if (!bRerunMD
|| !rerun_fr
.bStep
)
2740 /* increase the MD step number */
2745 cycles
= wallcycle_stop(wcycle
,ewcSTEP
);
2746 if (DOMAINDECOMP(cr
) && wcycle
)
2748 dd_cycles_add(cr
->dd
,cycles
,ddCyclStep
);
2751 if (step_rel
== wcycle_get_reset_counters(wcycle
) ||
2752 gs
.set
[eglsRESETCOUNTERS
] != 0)
2754 /* Reset all the counters related to performance over the run */
2755 reset_all_counters(fplog
,cr
,step
,&step_rel
,ir
,wcycle
,nrnb
,runtime
);
2756 wcycle_set_reset_counters(wcycle
,-1);
2757 /* Correct max_hours for the elapsed time */
2758 max_hours
-= run_time
/(60.0*60.0);
2759 bResetCountersHalfMaxH
= FALSE
;
2760 gs
.set
[eglsRESETCOUNTERS
] = 0;
2763 /* End of main MD loop */
2767 runtime_end(runtime
);
2769 if (bRerunMD
&& MASTER(cr
))
2774 if (!(cr
->duty
& DUTY_PME
))
2776 /* Tell the PME only node to finish */
2782 if (ir
->nstcalcenergy
> 0 && !bRerunMD
)
2784 print_ebin(outf
->fp_ene
,FALSE
,FALSE
,FALSE
,fplog
,step
,t
,
2785 eprAVER
,FALSE
,mdebin
,fcd
,groups
,&(ir
->opts
));
2793 if (ir
->nstlist
== -1 && nlh
.nns
> 0 && fplog
)
2795 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
)));
2796 fprintf(fplog
,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh
.ab
/nlh
.nns
);
2799 if (shellfc
&& fplog
)
2801 fprintf(fplog
,"Fraction of iterations that converged: %.2f %%\n",
2802 (nconverged
*100.0)/step_rel
);
2803 fprintf(fplog
,"Average number of force evaluations per MD step: %.2f\n\n",
2807 if (repl_ex_nst
> 0 && MASTER(cr
))
2809 print_replica_exchange_statistics(fplog
,repl_ex
);
2812 runtime
->nsteps_done
= step_rel
;