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
74 #include "mpelogging.h"
81 #include "compute_io.h"
83 #include "checkpoint.h"
84 #include "mtop_util.h"
98 /* The following two variables and the signal_handler function
99 * is used from pme.c as well
101 extern bool bGotTermSignal
, bGotUsr1Signal
;
103 static RETSIGTYPE
signal_handler(int n
)
107 bGotTermSignal
= TRUE
;
111 bGotUsr1Signal
= TRUE
;
118 gmx_integrator_t
*func
;
121 /* The array should match the eI array in include/types/enums.h */
122 const gmx_intp_t integrator
[eiNR
] = { {do_md
}, {do_steep
}, {do_cg
}, {do_md
}, {do_md
}, {do_nm
}, {do_lbfgs
}, {do_tpi
}, {do_tpi
}, {do_md
}, {do_md
}};
124 /* Static variables for temporary use with the deform option */
125 static int init_step_tpx
;
126 static matrix box_tpx
;
128 static tMPI_Thread_mutex_t box_mutex
=TMPI_THREAD_MUTEX_INITIALIZER
;
132 static void copy_coupling_state(t_state
*statea
,t_state
*stateb
,
133 gmx_ekindata_t
*ekinda
,gmx_ekindata_t
*ekindb
)
136 /* MRS note -- might be able to get rid of some of the arguments. Look over it when it's all debugged */
140 stateb
->natoms
= statea
->natoms
;
141 stateb
->ngtc
= statea
->ngtc
;
142 stateb
->veta
= statea
->veta
;
145 copy_mat(ekinda
->ekin
,ekindb
->ekin
);
146 for (i
=0; i
<stateb
->ngtc
; i
++)
148 ekindb
->tcstat
[i
].T
= ekinda
->tcstat
[i
].T
;
149 ekindb
->tcstat
[i
].Th
= ekinda
->tcstat
[i
].Th
;
150 copy_mat(ekinda
->tcstat
[i
].ekinh
,ekindb
->tcstat
[i
].ekinh
);
151 copy_mat(ekinda
->tcstat
[i
].ekinf
,ekindb
->tcstat
[i
].ekinf
);
152 ekindb
->tcstat
[i
].ekinscalef_nhc
= ekinda
->tcstat
[i
].ekinscalef_nhc
;
153 ekindb
->tcstat
[i
].ekinscaleh_nhc
= ekinda
->tcstat
[i
].ekinscaleh_nhc
;
154 ekindb
->tcstat
[i
].vscale_nhc
= ekinda
->tcstat
[i
].vscale_nhc
;
157 copy_rvecn(statea
->x
,stateb
->x
,0,stateb
->natoms
);
158 copy_rvecn(statea
->v
,stateb
->v
,0,stateb
->natoms
);
159 copy_mat(statea
->box
,stateb
->box
);
160 copy_mat(statea
->box_rel
,stateb
->box_rel
);
161 copy_mat(statea
->boxv
,stateb
->boxv
);
162 /* need extra term for the barostat */
163 ngtc_eff
= stateb
->ngtc
+1;
165 for (i
= 0; i
< ngtc_eff
; i
++)
167 for (j
=0; j
< NNHCHAIN
; j
++)
169 stateb
->nosehoover_xi
[i
*NNHCHAIN
+ j
] = statea
->nosehoover_xi
[i
*NNHCHAIN
+ j
];
170 stateb
->nosehoover_vxi
[i
*NNHCHAIN
+ j
] = statea
->nosehoover_vxi
[i
*NNHCHAIN
+ j
];
175 static void compute_globals(FILE *fplog
, gmx_global_stat_t gstat
, t_commrec
*cr
, t_inputrec
*ir
,
176 t_forcerec
*fr
, gmx_ekindata_t
*ekind
,
177 t_state
*state
, t_state
*state_global
, t_mdatoms
*mdatoms
,
178 t_nrnb
*nrnb
, t_vcm
*vcm
, gmx_wallcycle_t wcycle
,
179 gmx_enerdata_t
*enerd
,tensor force_vir
, tensor shake_vir
, tensor total_vir
,
180 tensor pres
, rvec mu_tot
, gmx_constr_t constr
,
181 real
*chkpt
,real
*terminate
, real
*terminate_now
,
182 int *nabnsb
, matrix box
, gmx_mtop_t
*top_global
, real
*pcurr
,
183 int natoms
, bool *bSumEkinhOld
, int flags
)
187 tensor corr_vir
,corr_pres
,shakeall_vir
;
188 bool bEner
,bPres
,bTemp
, bVV
;
189 bool bRerunMD
, bEkinAveVel
, bStopCM
, bGStat
, bNEMD
, bFirstHalf
, bIterate
,
190 bFirstIterate
, bCopyEkinh
, bReadEkin
,bScaleEkin
;
191 static real ekin
,temp
;
192 static real prescorr
,enercorr
,dvdlcorr
,vtemp
;
194 /* translate CGLO flags to booleans */
195 bRerunMD
= flags
& CGLO_RERUNMD
;
196 bEkinAveVel
= flags
& CGLO_EKINAVEVEL
;
197 bStopCM
= flags
& CGLO_STOPCM
;
198 bGStat
= flags
& CGLO_GSTAT
;
199 bNEMD
= flags
& CGLO_NEMD
;
200 bIterate
= flags
& CGLO_ITERATE
;
201 bFirstIterate
= flags
& CGLO_FIRSTITERATE
;
202 bCopyEkinh
= flags
& CGLO_COPYEKINH
;
203 bEner
= flags
& CGLO_ENERGY
;
204 bTemp
= flags
& CGLO_TEMPERATURE
;
205 bPres
= flags
& CGLO_PRESSURE
;
206 bReadEkin
= flags
& CGLO_READEKIN
;
207 bScaleEkin
= flags
& CGLO_SCALEEKIN
;
209 /* decide when to calculate temperature and pressure. */
211 /* in initalization, it sums the shake virial in vv, and to
212 sums ekinh_old in leapfrog (or if we are calculating ekinh_old for other reasons */
214 bVV
= (ir
->eI
== eiVV
);
216 /* ########## Kinetic energy ############## */
220 /* Non-equilibrium MD: this is parallellized, but only does communication
221 * when there really is NEMD.
224 if (PAR(cr
) && (bNEMD
))
226 accumulate_u(cr
,&(ir
->opts
),ekind
);
231 restore_ekinstate_from_state(cr
,ekind
,&state_global
->ekinstate
);
235 calc_ke_part(state
,&(ir
->opts
),mdatoms
,ekind
,nrnb
,bEkinAveVel
,bIterate
);
240 /* Calculate center of mass velocity if necessary, also parallellized */
241 if (bStopCM
&& !bRerunMD
&& bEner
)
243 calc_vcm_grp(fplog
,mdatoms
->start
,mdatoms
->homenr
,mdatoms
,
244 state
->x
,state
->v
,vcm
);
248 if (bTemp
|| bPres
| bEner
)
252 /* We will not sum ekinh_old,
253 * so signal that we still have to do it.
255 *bSumEkinhOld
= TRUE
;
261 GMX_MPE_LOG(ev_global_stat_start
);
262 global_stat(fplog
,gstat
,cr
,enerd
,force_vir
,shake_vir
,mu_tot
,
263 ir
,ekind
,constr
,vcm
,NULL
,NULL
,terminate
,top_global
,state
,
264 *bSumEkinhOld
,flags
);
265 GMX_MPE_LOG(ev_global_stat_finish
);
270 terminate_now
= terminate
;
273 wallcycle_stop(wcycle
,ewcMoveE
);
274 *bSumEkinhOld
= FALSE
;
278 if (!bNEMD
&& debug
&& bTemp
&& (vcm
->nr
> 0))
281 mdatoms
->start
,mdatoms
->start
+mdatoms
->homenr
,
282 state
->v
,vcm
->group_p
[0],
283 mdatoms
->massT
,mdatoms
->tmass
,ekind
->ekin
);
287 /* Do center of mass motion removal */
288 if (bStopCM
&& !bRerunMD
&& !bFirstHalf
) /* fix this! */
290 check_cm_grp(fplog
,vcm
,ir
,1);
291 do_stopcm_grp(fplog
,mdatoms
->start
,mdatoms
->homenr
,mdatoms
->cVCM
,
292 state
->x
,state
->v
,vcm
);
293 inc_nrnb(nrnb
,eNR_STOPCM
,mdatoms
->homenr
);
299 /* Sum the kinetic energies of the groups & calc temp */
300 enerd
->term
[F_TEMP
] = sum_ekin(&(ir
->opts
),ekind
,&(enerd
->term
[F_DKDL
]),
301 bEkinAveVel
||bReadEkin
,bIterate
,bScaleEkin
);
302 /* three main: VV with AveVel, vv with AveEkin, leap with AveEkin. Leap with AveVel is also
303 an option, but not supported now. Additionally, if we are doing iterations.
304 bCopyHalf: if TRUE, we simply copy the ekinh directly to ekin, multiplying be ekinscale_nhc.
305 bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale.
306 If FALSE, we average ekinh_old and ekinh*ekinscale_nhc.
307 bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
308 If FALSE, we go ahead and erase.
310 enerd
->term
[F_EKIN
] = trace(ekind
->ekin
);
313 /* ########## Long range energy information ###### */
317 calc_dispcorr(fplog
,ir
,fr
,0,top_global
,box
,state
->lambda
,
318 corr_pres
,corr_vir
,&prescorr
,&enercorr
,&dvdlcorr
);
321 if (bEner
&& bFirstIterate
)
323 enerd
->term
[F_DISPCORR
] = enercorr
;
324 enerd
->term
[F_EPOT
] += enercorr
;
325 enerd
->term
[F_DVDL
] += dvdlcorr
;
326 if (fr
->efep
!= efepNO
) {
327 enerd
->dvdl_lin
+= dvdlcorr
;
331 /* ########## Now pressure ############## */
335 m_add(force_vir
,shake_vir
,total_vir
);
337 /* Calculate pressure and apply LR correction if PPPM is used.
338 * Use the box from last timestep since we already called update().
341 enerd
->term
[F_PRES
] = calc_pres(fr
->ePBC
,ir
->nwall
,box
,ekind
->ekin
,total_vir
,pres
,
342 (fr
->eeltype
==eelPPPM
)?enerd
->term
[F_COUL_RECIP
]:0.0);
344 /* Calculate long range corrections to pressure and energy */
345 /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
346 and computes enerd->term[F_DISPCORR]. Also modifies the
347 total_vir and pres tesors */
349 m_add(total_vir
,corr_vir
,total_vir
);
350 m_add(pres
,corr_pres
,pres
);
351 enerd
->term
[F_PDISPCORR
] = prescorr
;
352 enerd
->term
[F_PRES
] += prescorr
;
353 *pcurr
= enerd
->term
[F_PRES
];
354 /* calculate temperature using virial */
355 enerd
->term
[F_VTEMP
] = calc_temp(trace(total_vir
),ir
->opts
.nrdf
[0]);
361 struct mdrunner_arglist
375 const char *dddlb_opt
;
388 int ret
; /* return value */
392 static void mdrunner_start_fn(void *arg
)
394 struct mdrunner_arglist
*mda
=(struct mdrunner_arglist
*)arg
;
395 struct mdrunner_arglist mc
=*mda
; /* copy the arg list to make sure
396 that it's thread-local. This doesn't
397 copy pointed-to items, of course,
398 but those are all const. */
399 t_commrec
*cr
; /* we need a local version of this */
401 t_filenm
*fnm
=dup_tfn(mc
.nfile
, mc
.fnm
);
403 cr
=init_par_threads(mc
.cr
);
410 mda
->ret
=mdrunner(fplog
, cr
, mc
.nfile
, mc
.fnm
, mc
.oenv
, mc
.bVerbose
,
411 mc
.bCompact
, mc
.nstglobalcomm
,
412 mc
.ddxyz
, mc
.dd_node_order
, mc
.rdd
,
413 mc
.rconstr
, mc
.dddlb_opt
, mc
.dlb_scale
,
414 mc
.ddcsx
, mc
.ddcsy
, mc
.ddcsz
, mc
.nstepout
, mc
.nmultisim
,
415 mc
.repl_ex_nst
, mc
.repl_ex_seed
, mc
.pforce
,
416 mc
.cpt_period
, mc
.max_hours
, mc
.Flags
);
421 int mdrunner_threads(int nthreads
,
422 FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
423 const output_env_t oenv
, bool bVerbose
,bool bCompact
,
425 ivec ddxyz
,int dd_node_order
,real rdd
,real rconstr
,
426 const char *dddlb_opt
,real dlb_scale
,
427 const char *ddcsx
,const char *ddcsy
,const char *ddcsz
,
428 int nstepout
,int nmultisim
,int repl_ex_nst
,
429 int repl_ex_seed
, real pforce
,real cpt_period
,
430 real max_hours
, unsigned long Flags
)
433 /* first check whether we even need to start tMPI */
436 ret
=mdrunner(fplog
, cr
, nfile
, fnm
, oenv
, bVerbose
, bCompact
,
438 ddxyz
, dd_node_order
, rdd
, rconstr
, dddlb_opt
, dlb_scale
,
439 ddcsx
, ddcsy
, ddcsz
, nstepout
, nmultisim
, repl_ex_nst
,
440 repl_ex_seed
, pforce
, cpt_period
, max_hours
, Flags
);
445 struct mdrunner_arglist mda
;
446 /* fill the data structure to pass as void pointer to thread start fn */
452 mda
.bVerbose
=bVerbose
;
453 mda
.bCompact
=bCompact
;
454 mda
.nstglobalcomm
=nstglobalcomm
;
455 mda
.ddxyz
[XX
]=ddxyz
[XX
];
456 mda
.ddxyz
[YY
]=ddxyz
[YY
];
457 mda
.ddxyz
[ZZ
]=ddxyz
[ZZ
];
458 mda
.dd_node_order
=dd_node_order
;
461 mda
.dddlb_opt
=dddlb_opt
;
462 mda
.dlb_scale
=dlb_scale
;
466 mda
.nstepout
=nstepout
;
467 mda
.nmultisim
=nmultisim
;
468 mda
.repl_ex_nst
=repl_ex_nst
;
469 mda
.repl_ex_seed
=repl_ex_seed
;
471 mda
.cpt_period
=cpt_period
;
472 mda
.max_hours
=max_hours
;
475 fprintf(stderr
, "Starting %d threads\n",nthreads
);
477 tMPI_Init_fn(nthreads
, mdrunner_start_fn
, (void*)(&mda
) );
481 gmx_comm("Multiple threads requested but not compiled with threads");
488 int mdrunner(FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
489 const output_env_t oenv
, bool bVerbose
,bool bCompact
,
491 ivec ddxyz
,int dd_node_order
,real rdd
,real rconstr
,
492 const char *dddlb_opt
,real dlb_scale
,
493 const char *ddcsx
,const char *ddcsy
,const char *ddcsz
,
494 int nstepout
,int nmultisim
,int repl_ex_nst
,int repl_ex_seed
,
495 real pforce
,real cpt_period
,real max_hours
,
498 double nodetime
=0,realtime
;
499 t_inputrec
*inputrec
;
506 gmx_mtop_t
*mtop
=NULL
;
507 t_mdatoms
*mdatoms
=NULL
;
511 gmx_pme_t
*pmedata
=NULL
;
512 gmx_vsite_t
*vsite
=NULL
;
514 int i
,m
,nChargePerturbed
=-1,status
,nalloc
;
516 gmx_wallcycle_t wcycle
;
517 bool bReadRNG
,bReadEkin
;
519 gmx_runtime_t runtime
;
521 gmx_large_int_t reset_counters
;
524 /* Essential dynamics */
525 if (opt2bSet("-ei",nfile
,fnm
))
527 /* Open input and output files, allocate space for ED data structure */
528 ed
= ed_open(nfile
,fnm
,cr
);
536 if (bVerbose
&& SIMMASTER(cr
))
538 fprintf(stderr
,"Getting Loaded...\n");
541 if (Flags
& MD_APPENDFILES
)
548 /* The master thread on the master node reads from disk,
549 * then distributes everything to the other processors.
552 list
= (SIMMASTER(cr
) && !(Flags
& MD_APPENDFILES
)) ? (LIST_SCALARS
| LIST_INPUTREC
) : 0;
555 init_parallel(fplog
, opt2fn_master("-s",nfile
,fnm
,cr
),cr
,
556 inputrec
,mtop
,state
,list
);
561 /* Read a file for a single processor */
563 init_single(fplog
,inputrec
,ftp2fn(efTPX
,nfile
,fnm
),mtop
,state
);
565 if (!EEL_PME(inputrec
->coulombtype
) || (Flags
& MD_PARTDEC
))
571 fcRegisterSteps(inputrec
->nsteps
,inputrec
->init_step
);
574 /* NMR restraints must be initialized before load_checkpoint,
575 * since with time averaging the history is added to t_state.
576 * For proper consistency check we therefore need to extend
578 * So the PME-only nodes (if present) will also initialize
579 * the distance restraints.
583 /* This needs to be called before read_checkpoint to extend the state */
584 init_disres(fplog
,mtop
,inputrec
,cr
,Flags
& MD_PARTDEC
,fcd
,state
);
586 if (gmx_mtop_ftype_count(mtop
,F_ORIRES
) > 0)
588 if (PAR(cr
) && !(Flags
& MD_PARTDEC
))
590 gmx_fatal(FARGS
,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
592 /* Orientation restraints */
595 init_orires(fplog
,mtop
,state
->x
,inputrec
,cr
->ms
,&(fcd
->orires
),
600 if (DEFORM(*inputrec
))
602 /* Store the deform reference box before reading the checkpoint */
605 copy_mat(state
->box
,box
);
609 gmx_bcast(sizeof(box
),box
,cr
);
611 /* Because we do not have the update struct available yet
612 * in which the reference values should be stored,
613 * we store them temporarily in static variables.
614 * This should be thread safe, since they are only written once
615 * and with identical values.
618 tMPI_Thread_mutex_lock(&box_mutex
);
620 init_step_tpx
= inputrec
->init_step
;
621 copy_mat(box
,box_tpx
);
623 tMPI_Thread_mutex_unlock(&box_mutex
);
627 if (opt2bSet("-cpi",nfile
,fnm
))
629 /* Check if checkpoint file exists before doing continuation.
630 * This way we can use identical input options for the first and subsequent runs...
632 if( gmx_fexist_master(opt2fn_master("-cpi",nfile
,fnm
,cr
),cr
) )
634 load_checkpoint(opt2fn_master("-cpi",nfile
,fnm
,cr
),fplog
,
635 cr
,Flags
& MD_PARTDEC
,ddxyz
,
636 inputrec
,state
,&bReadRNG
,&bReadEkin
,
637 (Flags
& MD_APPENDFILES
));
641 Flags
|= MD_READ_RNG
;
645 Flags
|= MD_READ_EKIN
;
650 if (MASTER(cr
) && (Flags
& MD_APPENDFILES
))
652 fplog
= gmx_log_open(ftp2fn(efLOG
,nfile
,fnm
),cr
,!(Flags
& MD_SEPPOT
),
658 copy_mat(state
->box
,box
);
663 gmx_bcast(sizeof(box
),box
,cr
);
666 if (bVerbose
&& SIMMASTER(cr
))
668 fprintf(stderr
,"Loaded with Money\n\n");
671 if (PAR(cr
) && !((Flags
& MD_PARTDEC
) || EI_TPI(inputrec
->eI
)))
673 cr
->dd
= init_domain_decomposition(fplog
,cr
,Flags
,ddxyz
,rdd
,rconstr
,
680 make_dd_communicators(fplog
,cr
,dd_node_order
);
682 /* Set overallocation to avoid frequent reallocation of arrays */
683 set_over_alloc_dd(TRUE
);
687 cr
->duty
= (DUTY_PP
| DUTY_PME
);
688 npme_major
= cr
->nnodes
;
690 if (inputrec
->ePBC
== epbcSCREW
)
693 "pbc=%s is only implemented with domain decomposition",
694 epbc_names
[inputrec
->ePBC
]);
700 /* After possible communicator splitting in make_dd_communicators.
701 * we can set up the intra/inter node communication.
703 gmx_setup_nodecomm(fplog
,cr
);
706 wcycle
= wallcycle_init(fplog
,cr
);
709 /* Master synchronizes its value of reset_counters with all nodes
710 * including PME only nodes */
711 reset_counters
= wcycle_get_reset_counters(wcycle
);
712 gmx_bcast_sim(sizeof(reset_counters
),&reset_counters
,cr
);
713 wcycle_set_reset_counters(wcycle
, reset_counters
);
718 if (cr
->duty
& DUTY_PP
)
720 /* For domain decomposition we allocate dynamically
721 * in dd_partition_system.
723 if (DOMAINDECOMP(cr
))
725 bcast_state_setup(cr
,state
);
735 bcast_state(cr
,state
,TRUE
);
739 /* Dihedral Restraints */
740 if (gmx_mtop_ftype_count(mtop
,F_DIHRES
) > 0)
742 init_dihres(fplog
,mtop
,inputrec
,fcd
);
745 /* Initiate forcerecord */
747 init_forcerec(fplog
,oenv
,fr
,fcd
,inputrec
,mtop
,cr
,box
,FALSE
,
748 opt2fn("-table",nfile
,fnm
),
749 opt2fn("-tablep",nfile
,fnm
),
750 opt2fn("-tableb",nfile
,fnm
),FALSE
,pforce
);
752 /* version for PCA_NOT_READ_NODE (see md.c) */
753 /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
754 "nofile","nofile","nofile",FALSE,pforce);
756 fr
->bSepDVDL
= ((Flags
& MD_SEPPOT
) == MD_SEPPOT
);
758 /* Initialize QM-MM */
761 init_QMMMrec(cr
,box
,mtop
,inputrec
,fr
);
764 /* Initialize the mdatoms structure.
765 * mdatoms is not filled with atom data,
766 * as this can not be done now with domain decomposition.
768 mdatoms
= init_mdatoms(fplog
,mtop
,inputrec
->efep
!=efepNO
);
770 /* Initialize the virtual site communication */
771 vsite
= init_vsite(mtop
,cr
);
773 calc_shifts(box
,fr
->shift_vec
);
775 /* With periodic molecules the charge groups should be whole at start up
776 * and the virtual sites should not be far from their proper positions.
778 if (!inputrec
->bContinuation
&& MASTER(cr
) &&
779 !(inputrec
->ePBC
!= epbcNONE
&& inputrec
->bPeriodicMols
))
781 /* Make molecules whole at start of run */
782 if (fr
->ePBC
!= epbcNONE
)
784 do_pbc_first_mtop(fplog
,inputrec
->ePBC
,box
,mtop
,state
->x
);
788 /* Correct initial vsite positions are required
789 * for the initial distribution in the domain decomposition
790 * and for the initial shell prediction.
792 construct_vsites_mtop(fplog
,vsite
,mtop
,state
->x
);
796 /* Initiate PPPM if necessary */
797 if (fr
->eeltype
== eelPPPM
)
799 if (mdatoms
->nChargePerturbed
)
801 gmx_fatal(FARGS
,"Free energy with %s is not implemented",
802 eel_names
[fr
->eeltype
]);
804 status
= gmx_pppm_init(fplog
,cr
,oenv
,FALSE
,TRUE
,box
,
805 getenv("GMXGHAT"),inputrec
, (Flags
& MD_REPRODUCIBLE
));
808 gmx_fatal(FARGS
,"Error %d initializing PPPM",status
);
812 if (EEL_PME(fr
->eeltype
))
814 ewaldcoeff
= fr
->ewaldcoeff
;
815 pmedata
= &fr
->pmedata
;
824 /* This is a PME only node */
826 /* We don't need the state */
829 ewaldcoeff
= calc_ewaldcoeff(inputrec
->rcoulomb
, inputrec
->ewald_rtol
);
833 /* Initiate PME if necessary,
834 * either on all nodes or on dedicated PME nodes only. */
835 if (EEL_PME(inputrec
->coulombtype
))
839 nChargePerturbed
= mdatoms
->nChargePerturbed
;
841 if (cr
->npmenodes
> 0)
843 /* The PME only nodes need to know nChargePerturbed */
844 gmx_bcast_sim(sizeof(nChargePerturbed
),&nChargePerturbed
,cr
);
846 if (cr
->duty
& DUTY_PME
)
848 status
= gmx_pme_init(pmedata
,cr
,npme_major
,inputrec
,
849 mtop
? mtop
->natoms
: 0,nChargePerturbed
,
850 (Flags
& MD_REPRODUCIBLE
));
853 gmx_fatal(FARGS
,"Error %d initializing PME",status
);
859 if (integrator
[inputrec
->eI
].func
== do_md
)
861 /* Turn on signal handling on all nodes */
863 * (A user signal from the PME nodes (if any)
864 * is communicated to the PP nodes.
866 if (getenv("GMX_NO_TERM") == NULL
)
870 fprintf(debug
,"Installing signal handler for SIGTERM\n");
872 signal(SIGTERM
,signal_handler
);
875 if (getenv("GMX_NO_USR1") == NULL
)
879 fprintf(debug
,"Installing signal handler for SIGUSR1\n");
881 signal(SIGUSR1
,signal_handler
);
886 if (cr
->duty
& DUTY_PP
)
888 if (inputrec
->ePull
!= epullNO
)
890 /* Initialize pull code */
891 init_pull(fplog
,inputrec
,nfile
,fnm
,mtop
,cr
,oenv
,
892 EI_DYNAMICS(inputrec
->eI
) && MASTER(cr
),Flags
);
895 constr
= init_constraints(fplog
,mtop
,inputrec
,ed
,state
,cr
);
897 if (DOMAINDECOMP(cr
))
899 dd_init_bondeds(fplog
,cr
->dd
,mtop
,vsite
,constr
,inputrec
,
900 Flags
& MD_DDBONDCHECK
,fr
->cginfo_mb
);
902 set_dd_parameters(fplog
,cr
->dd
,dlb_scale
,inputrec
,fr
,&ddbox
);
904 setup_dd_grid(fplog
,cr
->dd
);
907 /* Now do whatever the user wants us to do (how flexible...) */
908 integrator
[inputrec
->eI
].func(fplog
,cr
,nfile
,fnm
,
909 oenv
,bVerbose
,bCompact
,
912 nstepout
,inputrec
,mtop
,
914 mdatoms
,nrnb
,wcycle
,ed
,fr
,
915 repl_ex_nst
,repl_ex_seed
,
916 cpt_period
,max_hours
,
920 if (inputrec
->ePull
!= epullNO
)
922 finish_pull(fplog
,inputrec
->pull
);
928 gmx_pmeonly(*pmedata
,cr
,nrnb
,wcycle
,ewaldcoeff
,FALSE
,inputrec
);
931 if (EI_DYNAMICS(inputrec
->eI
) || EI_TPI(inputrec
->eI
))
933 /* Some timing stats */
936 if (runtime
.proc
== 0)
938 runtime
.proc
= runtime
.real
;
947 wallcycle_stop(wcycle
,ewcRUN
);
949 /* Finish up, write some stuff
950 * if rerunMD, don't write last frame again
952 finish_run(fplog
,cr
,ftp2fn(efSTO
,nfile
,fnm
),
953 inputrec
,nrnb
,wcycle
,&runtime
,
954 EI_DYNAMICS(inputrec
->eI
) && !MULTISIM(cr
));
956 /* Does what it says */
957 print_date_and_time(fplog
,cr
->nodeid
,"Finished mdrun",&runtime
);
959 /* Close logfile already here if we were appending to it */
960 if (MASTER(cr
) && (Flags
& MD_APPENDFILES
))
962 gmx_log_close(fplog
);
969 else if(bGotUsr1Signal
)
981 static void md_print_warning(const t_commrec
*cr
,FILE *fplog
,const char *buf
)
985 fprintf(stderr
,"\n%s\n",buf
);
989 fprintf(fplog
,"\n%s\n",buf
);
993 static bool done_iterating(const t_commrec
*cr
,FILE *fplog
, bool *bFirstIterate
, bool *bIterate
, real fom
, real
*newf
, int n
)
998 /* monitor convergence, and use a secant search to propose new
1001 The secant method computes x_{i+1} = x_{i} - f(x_{i}) * ---------------------
1002 f(x_{i}) - f(x_{i-1})
1004 The function we are trying to zero is fom-x, where fom is the
1005 "figure of merit" which is the pressure (or the veta value) we
1006 would get by putting in an old value of the pressure or veta into
1007 the incrementor function for the step or half step. I have
1008 verified that this gives the same answer as self consistent
1009 iteration, usually in many fewer steps, especially for small tau_p.
1011 We could possibly eliminate an iteration with proper use
1012 of the value from the previous step, but that would take a bit
1013 more bookkeeping, especially for veta, since tests indicate the
1014 function of veta on the last step is not sufficiently close to
1015 guarantee convergence this step. This is
1016 good enough for now. On my tests, I could use tau_p down to
1017 0.02, which is smaller that would ever be necessary in
1018 practice. Generally, 3-5 iterations will be sufficient */
1020 /* Definitions for convergence. Could be optimized . . . */
1022 #define CONVERGEITER 0.000000001
1024 #define CONVERGEITER 0.0001
1026 #define MAXITERCONST 200
1028 /* used to escape out of cyclic traps because of limited numberical precision */
1030 #define FALLBACK 1000
1032 static real f
,fprev
,x
,xprev
;
1034 static real allrelerr
[MAXITERCONST
+2];
1042 *bFirstIterate
= FALSE
;
1050 f
= fom
-x
; /* we want to zero this difference */
1051 if ((iter_i
> 1) && (iter_i
< MAXITERCONST
))
1059 *newf
= x
- (x
-xprev
)*(f
)/(f
-fprev
);
1064 /* just use self-consistent iteration the first step to initialize, or
1065 if it's not converging (which happens occasionally -- need to investigate why) */
1069 /* Consider a slight shortcut allowing us to exit one sooner -- we check the
1070 difference between the closest of x and xprev to the new
1071 value. To be 100% certain, we should check the difference between
1072 the last result, and the previous result, or
1074 relerr = (fabs((x-xprev)/fom));
1076 but this is pretty much never necessary under typical conditions.
1077 Checking numerically, it seems to lead to almost exactly the same
1078 trajectories, but there are small differences out a few decimal
1079 places in the pressure, and eventually in the v_eta, but it could
1082 if (fabs(*newf-x) < fabs(*newf - xprev)) { xmin = x;} else { xmin = xprev;}
1083 relerr = (fabs((*newf-xmin) / *newf));
1086 relerr
= (fabs((f
-fprev
)/fom
));
1088 allrelerr
[iter_i
] = relerr
;
1094 fprintf(debug
,"Iterating NPT constraints #%i: %6i %20.12f%14.6g%20.12f\n",n
,iter_i
,fom
,relerr
,*newf
);
1097 if ((relerr
< CONVERGEITER
) || (fom
==0))
1102 fprintf(debug
,"Iterating NPT constraints #%i: CONVERGED\n",n
);
1106 if (iter_i
> MAXITERCONST
)
1109 for (i
=0;i
<CYCLEMAX
;i
++) {
1110 if (allrelerr
[(MAXITERCONST
-2*CYCLEMAX
)-i
] == allrelerr
[iter_i
-1]) {
1112 if (relerr
> CONVERGEITER
*FALLBACK
) {incycle
= FALSE
; break;}
1117 /* we are trapped in a numerical attractor, and can't converge any more, and are close to the final result.
1118 Better to give up here and proceed with a warning than have the simulation die.
1120 sprintf(buf
,"numerical convergence issues with NPT, relative error only %10.5g, continuing\n",relerr
);
1121 md_print_warning(cr
,fplog
,buf
);
1124 gmx_fatal(FARGS
,"Could not converge NPT constraints\n");
1137 static void check_nst_param(FILE *fplog
,t_commrec
*cr
,
1138 const char *desc_nst
,int nst
,
1139 const char *desc_p
,int *p
)
1143 if (*p
> 0 && *p
% nst
!= 0)
1145 /* Round up to the next multiple of nst */
1146 *p
= ((*p
)/nst
+ 1)*nst
;
1147 sprintf(buf
,"NOTE: %s changes %s to %d\n",desc_nst
,desc_p
,*p
);
1148 md_print_warning(cr
,fplog
,buf
);
1152 static void reset_all_counters(FILE *fplog
,t_commrec
*cr
,
1153 gmx_large_int_t step
,
1154 gmx_large_int_t
*step_rel
,t_inputrec
*ir
,
1155 gmx_wallcycle_t wcycle
,t_nrnb
*nrnb
,
1156 gmx_runtime_t
*runtime
)
1158 char buf
[STRLEN
],sbuf
[22];
1160 /* Reset all the counters related to performance over the run */
1161 sprintf(buf
,"Step %s: resetting all time and cycle counters\n",
1162 gmx_step_str(step
,sbuf
));
1163 md_print_warning(cr
,fplog
,buf
);
1165 wallcycle_stop(wcycle
,ewcRUN
);
1166 wallcycle_reset_all(wcycle
);
1167 if (DOMAINDECOMP(cr
))
1169 reset_dd_statistics_counters(cr
->dd
);
1172 ir
->init_step
+= *step_rel
;
1173 ir
->nsteps
-= *step_rel
;
1175 wallcycle_start(wcycle
,ewcRUN
);
1176 runtime_start(runtime
);
1177 print_date_and_time(fplog
,cr
->nodeid
,"Restarted time",runtime
);
1180 static int check_nstglobalcomm(FILE *fplog
,t_commrec
*cr
,
1181 int nstglobalcomm
,t_inputrec
*ir
)
1185 if (!EI_DYNAMICS(ir
->eI
))
1190 if (nstglobalcomm
== -1)
1192 if (ir
->nstcalcenergy
== 0 && ir
->nstlist
== 0)
1195 if (ir
->nstenergy
> 0 && ir
->nstenergy
< nstglobalcomm
)
1197 nstglobalcomm
= ir
->nstenergy
;
1202 /* We assume that if nstcalcenergy > nstlist,
1203 * nstcalcenergy is a multiple of nstlist.
1205 if (ir
->nstcalcenergy
== 0 ||
1206 (ir
->nstlist
> 0 && ir
->nstlist
< ir
->nstcalcenergy
))
1208 nstglobalcomm
= ir
->nstlist
;
1212 nstglobalcomm
= ir
->nstcalcenergy
;
1218 if (ir
->nstlist
> 0 &&
1219 nstglobalcomm
> ir
->nstlist
&& nstglobalcomm
% ir
->nstlist
!= 0)
1221 nstglobalcomm
= (nstglobalcomm
/ ir
->nstlist
)*ir
->nstlist
;
1222 sprintf(buf
,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm
);
1223 md_print_warning(cr
,fplog
,buf
);
1225 if (nstglobalcomm
> ir
->nstcalcenergy
)
1227 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
1228 "nstcalcenergy",&ir
->nstcalcenergy
);
1231 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
1232 "nstenergy",&ir
->nstenergy
);
1234 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
1235 "nstlog",&ir
->nstlog
);
1238 if (ir
->comm_mode
!= ecmNO
&& ir
->nstcomm
< nstglobalcomm
)
1240 sprintf(buf
,"WARNING: Changing nstcomm from %d to %d\n",
1241 ir
->nstcomm
,nstglobalcomm
);
1242 md_print_warning(cr
,fplog
,buf
);
1243 ir
->nstcomm
= nstglobalcomm
;
1246 return nstglobalcomm
;
1249 static void check_ir_old_tpx_versions(t_commrec
*cr
,FILE *fplog
,
1250 t_inputrec
*ir
,gmx_mtop_t
*mtop
)
1252 /* Check required for old tpx files */
1253 if (IR_TWINRANGE(*ir
) && ir
->nstlist
> 1 &&
1254 ir
->nstcalcenergy
% ir
->nstlist
!= 0)
1256 md_print_warning(cr
,fplog
,"Old tpr file with twin-range settings: modifiying energy calculation and/or T/P-coupling frequencies");
1258 if (gmx_mtop_ftype_count(mtop
,F_CONSTR
) +
1259 gmx_mtop_ftype_count(mtop
,F_CONSTRNC
) > 0 &&
1260 ir
->eConstrAlg
== econtSHAKE
)
1262 md_print_warning(cr
,fplog
,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect");
1263 if (ir
->epc
!= epcNO
)
1265 gmx_fatal(FARGS
,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
1268 check_nst_param(fplog
,cr
,"nstlist",ir
->nstlist
,
1269 "nstcalcenergy",&ir
->nstcalcenergy
);
1271 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
1272 "nstenergy",&ir
->nstenergy
);
1273 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
1274 "nstlog",&ir
->nstlog
);
1275 if (ir
->efep
!= efepNO
)
1277 check_nst_param(fplog
,cr
,"nstdhdl",ir
->nstdhdl
,
1278 "nstenergy",&ir
->nstenergy
);
1283 bool bGStatEveryStep
;
1284 gmx_large_int_t step_ns
;
1285 gmx_large_int_t step_nscheck
;
1286 gmx_large_int_t nns
;
1296 static void reset_nlistheuristics(gmx_nlheur_t
*nlh
,gmx_large_int_t step
)
1300 nlh
->step_nscheck
= step
;
1303 static void init_nlistheuristics(gmx_nlheur_t
*nlh
,
1304 bool bGStatEveryStep
,gmx_large_int_t step
)
1306 nlh
->bGStatEveryStep
= bGStatEveryStep
;
1313 reset_nlistheuristics(nlh
,step
);
1316 static void update_nliststatistics(gmx_nlheur_t
*nlh
,gmx_large_int_t step
)
1318 gmx_large_int_t nl_lt
;
1319 char sbuf
[22],sbuf2
[22];
1321 /* Determine the neighbor list life time */
1322 nl_lt
= step
- nlh
->step_ns
;
1325 fprintf(debug
,"%d atoms beyond ns buffer, updating neighbor list after %s steps\n",nlh
->nabnsb
,gmx_step_str(nl_lt
,sbuf
));
1329 nlh
->s2
+= nl_lt
*nl_lt
;
1330 nlh
->ab
+= nlh
->nabnsb
;
1331 if (nlh
->lt_runav
== 0)
1333 nlh
->lt_runav
= nl_lt
;
1334 /* Initialize the fluctuation average
1335 * such that at startup we check after 0 steps.
1337 nlh
->lt_runav2
= sqr(nl_lt
/2.0);
1339 /* Running average with 0.9 gives an exp. history of 9.5 */
1340 nlh
->lt_runav2
= 0.9*nlh
->lt_runav2
+ 0.1*sqr(nlh
->lt_runav
- nl_lt
);
1341 nlh
->lt_runav
= 0.9*nlh
->lt_runav
+ 0.1*nl_lt
;
1342 if (nlh
->bGStatEveryStep
)
1344 /* Always check the nlist validity */
1345 nlh
->step_nscheck
= step
;
1349 /* We check after: <life time> - 2*sigma
1350 * The factor 2 is quite conservative,
1351 * but we assume that with nstlist=-1 the user
1352 * prefers exact integration over performance.
1354 nlh
->step_nscheck
= step
1355 + (int)(nlh
->lt_runav
- 2.0*sqrt(nlh
->lt_runav2
)) - 1;
1359 fprintf(debug
,"nlist life time %s run av. %4.1f sig %3.1f check %s check with -gcom %d\n",
1360 gmx_step_str(nl_lt
,sbuf
),nlh
->lt_runav
,sqrt(nlh
->lt_runav2
),
1361 gmx_step_str(nlh
->step_nscheck
-step
+1,sbuf2
),
1362 (int)(nlh
->lt_runav
- 2.0*sqrt(nlh
->lt_runav2
)));
1366 static void set_nlistheuristics(gmx_nlheur_t
*nlh
,bool bReset
,gmx_large_int_t step
)
1372 reset_nlistheuristics(nlh
,step
);
1376 update_nliststatistics(nlh
,step
);
1379 nlh
->step_ns
= step
;
1380 /* Initialize the cumulative coordinate scaling matrix */
1381 clear_mat(nlh
->scale_tot
);
1382 for(d
=0; d
<DIM
; d
++)
1384 nlh
->scale_tot
[d
][d
] = 1.0;
1388 double do_md(FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
1389 const output_env_t oenv
, bool bVerbose
,bool bCompact
,
1391 gmx_vsite_t
*vsite
,gmx_constr_t constr
,
1392 int stepout
,t_inputrec
*ir
,
1393 gmx_mtop_t
*top_global
,
1395 t_state
*state_global
,
1397 t_nrnb
*nrnb
,gmx_wallcycle_t wcycle
,
1398 gmx_edsam_t ed
,t_forcerec
*fr
,
1399 int repl_ex_nst
,int repl_ex_seed
,
1400 real cpt_period
,real max_hours
,
1401 unsigned long Flags
,
1402 gmx_runtime_t
*runtime
)
1404 int fp_trn
=0,fp_xtc
=0;
1405 ener_file_t fp_ene
=NULL
;
1406 gmx_large_int_t step
,step_rel
;
1408 FILE *fp_dhdl
=NULL
,*fp_field
=NULL
;
1411 bool bGStatEveryStep
,bGStat
,bNstEner
,bCalcPres
,bCalcEner
;
1412 bool bNS
,bNStList
,bSimAnn
,bStopCM
,bRerunMD
,bNotLastFrame
=FALSE
,
1413 bFirstStep
,bStateFromTPX
,bInitStep
,bLastStep
,bEkinAveVel
,
1414 bBornRadii
,bStartingFromCpt
,bVVAveVel
,bVVAveEkin
;
1416 bool bNEMD
,do_ene
,do_log
,do_verbose
,bRerunWarnNoV
=TRUE
,
1417 bForceUpdate
=FALSE
,bX
,bV
,bF
,bXTC
,bCPT
;
1419 int force_flags
,cglo_flags
;
1420 tensor force_vir
,shake_vir
,total_vir
,tmp_vir
,pres
;
1425 matrix
*scale_tot
,pcoupl_mu
,M
,ebox
;
1427 t_trxframe rerun_fr
;
1428 gmx_repl_ex_t repl_ex
=NULL
;
1430 /* Booleans (disguised as a reals) to checkpoint and terminate mdrun */
1431 real chkpt
=0,terminate
=0,terminate_now
=0;
1433 gmx_localtop_t
*top
;
1434 t_mdebin
*mdebin
=NULL
;
1435 t_state
*state
=NULL
;
1436 rvec
*f_global
=NULL
;
1439 gmx_enerdata_t
*enerd
;
1441 gmx_global_stat_t gstat
;
1442 gmx_update_t upd
=NULL
;
1443 t_graph
*graph
=NULL
;
1446 gmx_groups_t
*groups
;
1447 gmx_ekindata_t
*ekind
, *ekind_save
;
1448 gmx_shellfc_t shellfc
;
1449 int count
,nconverged
=0;
1453 bool bTCR
=FALSE
,bConverged
=TRUE
,bOK
,bSumEkinhOld
,bExchanged
;
1455 bool bVV
,bIterate
,bIterateFirstHalf
,bFirstIterate
,bTemp
,bPres
,bTrotter
;
1456 real temp0
,mu_aver
=0,dvdl
;
1458 atom_id
*grpindex
=NULL
;
1460 t_coupl_rec
*tcr
=NULL
;
1461 rvec
*xcopy
=NULL
,*vcopy
=NULL
,*cbuf
=NULL
;
1462 matrix boxcopy
,lastbox
;
1464 real fom
,oldfom
,veta_save
,pcurr
,scalevir
,tracevir
;
1467 int reset_counters
=-1;
1468 real last_conserved
= 0;
1473 char sbuf
[22],sbuf2
[22];
1474 bool bHandledSignal
=FALSE
;
1476 /* Temporary addition for FAHCORE checkpointing */
1480 /* Check for special mdrun options */
1481 bRerunMD
= (Flags
& MD_RERUN
);
1482 bIonize
= (Flags
& MD_IONIZE
);
1483 bFFscan
= (Flags
& MD_FFSCAN
);
1484 bAppend
= (Flags
& MD_APPENDFILES
);
1485 /* md-vv uses averaged half step velocities to determine temperature unless defined otherwise by GMX_EKIN_AVE_EKIN;
1486 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
1487 bVV
= (ir
->eI
==eiVV
);
1489 bEkinAveVel
= (getenv("GMX_EKIN_AVE_EKIN")==NULL
);
1491 bEkinAveVel
= (getenv("GMX_EKIN_AVE_VEL")!=NULL
);
1493 bVVAveVel
= bVV
&& bEkinAveVel
;
1494 bVVAveEkin
= bVV
&& !bEkinAveVel
;
1496 if (bVV
) /* to store the initial velocities while computing virial */
1498 snew(cbuf
,top_global
->natoms
);
1500 /* all the iteratative cases - only if there are constraints */
1501 bIterate
= ((IR_NPT_TROTTER(ir
)) && (constr
) && (!bRerunMD
));
1502 bTrotter
= (bVV
&& (IR_NPT_TROTTER(ir
) || (IR_NVT_TROTTER(ir
))));
1506 /* Since we don't know if the frames read are related in any way,
1507 * rebuild the neighborlist at every step.
1510 ir
->nstcalcenergy
= 1;
1514 check_ir_old_tpx_versions(cr
,fplog
,ir
,top_global
);
1516 nstglobalcomm
= check_nstglobalcomm(fplog
,cr
,nstglobalcomm
,ir
);
1517 bGStatEveryStep
= (nstglobalcomm
== 1);
1519 if (!bGStatEveryStep
&& ir
->nstlist
== -1 && fplog
!= NULL
)
1522 "To reduce the energy communication with nstlist = -1\n"
1523 "the neighbor list validity should not be checked at every step,\n"
1524 "this means that exact integration is not guaranteed.\n"
1525 "The neighbor list validity is checked after:\n"
1526 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
1527 "In most cases this will result in exact integration.\n"
1528 "This reduces the energy communication by a factor of 2 to 3.\n"
1529 "If you want less energy communication, set nstlist > 3.\n\n");
1532 if (bRerunMD
|| bFFscan
)
1536 groups
= &top_global
->groups
;
1538 /* Initial values */
1539 init_md(fplog
,cr
,ir
,oenv
,&t
,&t0
,&state_global
->lambda
,&lam0
,
1540 nrnb
,top_global
,&upd
,
1541 nfile
,fnm
,&fp_trn
,&fp_xtc
,&fp_ene
,&fn_cpt
,
1542 &fp_dhdl
,&fp_field
,&mdebin
,
1543 force_vir
,shake_vir
,mu_tot
,&bNEMD
,&bSimAnn
,&vcm
,Flags
);
1545 clear_mat(total_vir
);
1547 /* Energy terms and groups */
1549 init_enerdata(top_global
->groups
.grps
[egcENER
].nr
,ir
->n_flambda
,enerd
);
1550 if (DOMAINDECOMP(cr
))
1556 snew(f
,top_global
->natoms
);
1559 /* Kinetic energy data */
1561 init_ekindata(fplog
,top_global
,&(ir
->opts
),ekind
);
1562 /* needed for iteration of constraints */
1564 init_ekindata(fplog
,top_global
,&(ir
->opts
),ekind_save
);
1565 /* Copy the cos acceleration to the groups struct */
1566 ekind
->cosacc
.cos_accel
= ir
->cos_accel
;
1568 gstat
= global_stat_init(ir
);
1571 /* Check for polarizable models and flexible constraints */
1572 shellfc
= init_shell_flexcon(fplog
,
1573 top_global
,n_flexible_constraints(constr
),
1574 (ir
->bContinuation
||
1575 (DOMAINDECOMP(cr
) && !MASTER(cr
))) ?
1576 NULL
: state_global
->x
);
1581 tMPI_Thread_mutex_lock(&box_mutex
);
1583 set_deform_reference_box(upd
,init_step_tpx
,box_tpx
);
1585 tMPI_Thread_mutex_unlock(&box_mutex
);
1590 double io
= compute_io(ir
,top_global
->natoms
,groups
,mdebin
->ebin
->nener
,1);
1591 if ((io
> 2000) && MASTER(cr
))
1593 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
1597 if (DOMAINDECOMP(cr
)) {
1598 top
= dd_init_local_top(top_global
);
1601 dd_init_local_state(cr
->dd
,state_global
,state
);
1603 if (DDMASTER(cr
->dd
) && ir
->nstfout
) {
1604 snew(f_global
,state_global
->natoms
);
1608 /* Initialize the particle decomposition and split the topology */
1609 top
= split_system(fplog
,top_global
,ir
,cr
);
1611 pd_cg_range(cr
,&fr
->cg0
,&fr
->hcg
);
1612 pd_at_range(cr
,&a0
,&a1
);
1614 top
= gmx_mtop_generate_local_top(top_global
,ir
);
1617 a1
= top_global
->natoms
;
1620 state
= partdec_init_local_state(cr
,state_global
);
1623 atoms2md(top_global
,ir
,0,NULL
,a0
,a1
-a0
,mdatoms
);
1626 set_vsite_top(vsite
,top
,mdatoms
,cr
);
1629 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
) {
1630 graph
= mk_graph(fplog
,&(top
->idef
),0,top_global
->natoms
,FALSE
,FALSE
);
1634 make_local_shells(cr
,mdatoms
,shellfc
);
1637 if (ir
->pull
&& PAR(cr
)) {
1638 dd_make_local_pull_groups(NULL
,ir
->pull
,mdatoms
);
1642 if (DOMAINDECOMP(cr
))
1644 /* Distribute the charge groups over the nodes from the master node */
1645 dd_partition_system(fplog
,ir
->init_step
,cr
,TRUE
,1,
1646 state_global
,top_global
,ir
,
1647 state
,&f
,mdatoms
,top
,fr
,
1648 vsite
,shellfc
,constr
,
1652 /* If not DD, copy gb data */
1653 if(ir
->implicit_solvent
&& !DOMAINDECOMP(cr
))
1655 make_local_gb(cr
,fr
->born
,ir
->gb_algorithm
);
1658 update_mdatoms(mdatoms
,state
->lambda
);
1662 /* Update mdebin with energy history if appending to output files */
1663 if ( Flags
& MD_APPENDFILES
)
1665 restore_energyhistory_from_state(mdebin
,&state_global
->enerhist
);
1667 /* Set the initial energy history in state to zero by updating once */
1668 update_energyhistory(&state_global
->enerhist
,mdebin
);
1671 if ((state
->flags
& (1<<estLD_RNG
)) && (Flags
& MD_READ_RNG
)) {
1672 /* Set the random state if we read a checkpoint file */
1673 set_stochd_state(upd
,state
);
1676 /* Initialize constraints */
1678 if (!DOMAINDECOMP(cr
))
1679 set_constraints(constr
,top
,ir
,mdatoms
,cr
);
1682 /* Check whether we have to GCT stuff */
1683 bTCR
= ftp2bSet(efGCT
,nfile
,fnm
);
1686 fprintf(stderr
,"Will do General Coupling Theory!\n");
1688 gnx
= top_global
->mols
.nr
;
1690 for(i
=0; (i
<gnx
); i
++) {
1695 if (repl_ex_nst
> 0 && MASTER(cr
))
1696 repl_ex
= init_replica_exchange(fplog
,cr
->ms
,state_global
,ir
,
1697 repl_ex_nst
,repl_ex_seed
);
1699 if (!ir
->bContinuation
&& !bRerunMD
)
1701 if (mdatoms
->cFREEZE
&& (state
->flags
& (1<<estV
)))
1703 /* Set the velocities of frozen particles to zero */
1704 for(i
=mdatoms
->start
; i
<mdatoms
->start
+mdatoms
->homenr
; i
++)
1706 for(m
=0; m
<DIM
; m
++)
1708 if (ir
->opts
.nFreeze
[mdatoms
->cFREEZE
[i
]][m
])
1718 /* Constrain the initial coordinates and velocities */
1719 do_constrain_first(fplog
,constr
,ir
,mdatoms
,state
,f
,
1720 graph
,cr
,nrnb
,fr
,top
,shake_vir
);
1724 /* Construct the virtual sites for the initial configuration */
1725 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,NULL
,
1726 top
->idef
.iparams
,top
->idef
.il
,
1727 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
1733 /* would need to modify if vv starts with full time step */
1734 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
1735 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
1736 constr
,&chkpt
,&terminate
,&terminate_now
,&(nlh
.nabnsb
),state
->box
,
1737 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
1739 (bVV
? CGLO_PRESSURE
:0) |
1740 (CGLO_TEMPERATURE
) |
1741 (bRerunMD
? CGLO_RERUNMD
:0) |
1742 (bEkinAveVel
? CGLO_EKINAVEVEL
:0) |
1744 ((Flags
& MD_READ_EKIN
) ? CGLO_READEKIN
:0)));
1745 /* I'm assuming we need global communication the first time */
1746 /* Calculate the initial half step temperature, and save the ekinh_old */
1747 if (!(Flags
& MD_STARTFROMCPT
))
1749 for(i
=0; (i
<ir
->opts
.ngtc
); i
++)
1751 copy_mat(ekind
->tcstat
[i
].ekinh
,ekind
->tcstat
[i
].ekinh_old
);
1754 temp0
= enerd
->term
[F_TEMP
];
1756 /* Initiate data for the special cases */
1759 bufstate
= init_bufstate(state
->natoms
,(ir
->opts
.ngtc
+1)*NNHCHAIN
); /* extra state for barostat */
1764 snew(xcopy
,state
->natoms
);
1765 snew(vcopy
,state
->natoms
);
1766 copy_rvecn(state
->x
,xcopy
,0,state
->natoms
);
1767 copy_rvecn(state
->v
,vcopy
,0,state
->natoms
);
1768 copy_mat(state
->box
,boxcopy
);
1771 /* need to make an initial call to get the Trotter variables set. */
1772 trotter_seq
= init_trotter(ir
,state
,&MassQ
,bTrotter
);
1776 if (constr
&& !ir
->bContinuation
&& ir
->eConstrAlg
== econtLINCS
)
1779 "RMS relative constraint deviation after constraining: %.2e\n",
1780 constr_rmsd(constr
,FALSE
));
1782 fprintf(fplog
,"Initial temperature: %g K\n",enerd
->term
[F_TEMP
]);
1785 fprintf(stderr
,"starting md rerun '%s', reading coordinates from"
1786 " input trajectory '%s'\n\n",
1787 *(top_global
->name
),opt2fn("-rerun",nfile
,fnm
));
1790 fprintf(stderr
,"Calculated time to finish depends on nsteps from "
1791 "run input file,\nwhich may not correspond to the time "
1792 "needed to process input trajectory.\n\n");
1797 fprintf(stderr
,"starting mdrun '%s'\n",
1798 *(top_global
->name
));
1799 if (ir
->init_step
> 0)
1801 fprintf(stderr
,"%s steps, %8.1f ps (continuing from step %s, %8.1f ps).\n",
1802 gmx_step_str(ir
->init_step
+ir
->nsteps
,sbuf
),
1803 (ir
->init_step
+ir
->nsteps
)*ir
->delta_t
,
1804 gmx_step_str(ir
->init_step
,sbuf2
),
1805 ir
->init_step
*ir
->delta_t
);
1809 fprintf(stderr
,"%s steps, %8.1f ps.\n",
1810 gmx_step_str(ir
->nsteps
,sbuf
),ir
->nsteps
*ir
->delta_t
);
1813 fprintf(fplog
,"\n");
1816 /* Set and write start time */
1817 runtime_start(runtime
);
1818 print_date_and_time(fplog
,cr
->nodeid
,"Started mdrun",runtime
);
1819 wallcycle_start(wcycle
,ewcRUN
);
1821 fprintf(fplog
,"\n");
1823 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
1825 chkpt_ret
=fcCheckPointParallel( cr
->nodeid
,
1827 if ( chkpt_ret
== 0 )
1828 gmx_fatal( 3,__FILE__
,__LINE__
, "Checkpoint error on step %d\n", 0 );
1832 /***********************************************************
1834 * Loop over MD steps
1836 ************************************************************/
1838 /* if rerunMD then read coordinates and velocities from input trajectory */
1841 if (getenv("GMX_FORCE_UPDATE"))
1843 bForceUpdate
= TRUE
;
1846 bNotLastFrame
= read_first_frame(oenv
,&status
,
1847 opt2fn("-rerun",nfile
,fnm
),
1848 &rerun_fr
,TRX_NEED_X
| TRX_READ_V
);
1849 if (rerun_fr
.natoms
!= top_global
->natoms
)
1852 "Number of atoms in trajectory (%d) does not match the "
1853 "run input file (%d)\n",
1854 rerun_fr
.natoms
,top_global
->natoms
);
1856 if (ir
->ePBC
!= epbcNONE
)
1860 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
);
1862 if (max_cutoff2(ir
->ePBC
,rerun_fr
.box
) < sqr(fr
->rlistlong
))
1864 gmx_fatal(FARGS
,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr
.step
,rerun_fr
.time
);
1867 /* Set the shift vectors.
1868 * Necessary here when have a static box different from the tpr box.
1870 calc_shifts(rerun_fr
.box
,fr
->shift_vec
);
1874 /* loop over MD steps or if rerunMD to end of input trajectory */
1876 /* Skip the first Nose-Hoover integration when we get the state from tpx */
1877 bStateFromTPX
= !opt2bSet("-cpi",nfile
,fnm
);
1878 bInitStep
= bFirstStep
&& (bStateFromTPX
|| bVV
);
1879 bStartingFromCpt
= (Flags
& MD_STARTFROMCPT
) && bInitStep
;
1881 bSumEkinhOld
= FALSE
;
1884 step
= ir
->init_step
;
1887 if (ir
->nstlist
== -1)
1889 init_nlistheuristics(&nlh
,bGStatEveryStep
,step
);
1892 bLastStep
= (bRerunMD
|| step_rel
> ir
->nsteps
);
1893 while (!bLastStep
|| (bRerunMD
&& bNotLastFrame
)) {
1895 wallcycle_start(wcycle
,ewcSTEP
);
1897 GMX_MPE_LOG(ev_timestep1
);
1900 if (rerun_fr
.bStep
) {
1901 step
= rerun_fr
.step
;
1902 step_rel
= step
- ir
->init_step
;
1904 if (rerun_fr
.bTime
) {
1914 bLastStep
= (step_rel
== ir
->nsteps
);
1915 t
= t0
+ step
*ir
->delta_t
;
1918 if (ir
->efep
!= efepNO
)
1920 if (bRerunMD
&& rerun_fr
.bLambda
&& (ir
->delta_lambda
!=0))
1922 state_global
->lambda
= rerun_fr
.lambda
;
1926 state_global
->lambda
= lam0
+ step
*ir
->delta_lambda
;
1928 state
->lambda
= state_global
->lambda
;
1929 bDoDHDL
= do_per_step(step
,ir
->nstdhdl
);
1934 update_annealing_target_temp(&(ir
->opts
),t
);
1939 if (!(DOMAINDECOMP(cr
) && !MASTER(cr
)))
1941 for(i
=0; i
<state_global
->natoms
; i
++)
1943 copy_rvec(rerun_fr
.x
[i
],state_global
->x
[i
]);
1947 for(i
=0; i
<state_global
->natoms
; i
++)
1949 copy_rvec(rerun_fr
.v
[i
],state_global
->v
[i
]);
1954 for(i
=0; i
<state_global
->natoms
; i
++)
1956 clear_rvec(state_global
->v
[i
]);
1960 fprintf(stderr
,"\nWARNING: Some frames do not contain velocities.\n"
1961 " Ekin, temperature and pressure are incorrect,\n"
1962 " the virial will be incorrect when constraints are present.\n"
1964 bRerunWarnNoV
= FALSE
;
1968 copy_mat(rerun_fr
.box
,state_global
->box
);
1969 copy_mat(state_global
->box
,state
->box
);
1971 if (vsite
&& (Flags
& MD_RERUN_VSITE
))
1973 if (DOMAINDECOMP(cr
))
1975 gmx_fatal(FARGS
,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
1979 /* Following is necessary because the graph may get out of sync
1980 * with the coordinates if we only have every N'th coordinate set
1982 mk_mshift(fplog
,graph
,fr
->ePBC
,state
->box
,state
->x
);
1983 shift_self(graph
,state
->box
,state
->x
);
1985 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,state
->v
,
1986 top
->idef
.iparams
,top
->idef
.il
,
1987 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
1990 unshift_self(graph
,state
->box
,state
->x
);
1995 /* Stop Center of Mass motion */
1996 bStopCM
= (ir
->comm_mode
!= ecmNO
&& do_per_step(step
,ir
->nstcomm
));
1998 /* Copy back starting coordinates in case we're doing a forcefield scan */
2001 for(ii
=0; (ii
<state
->natoms
); ii
++)
2003 copy_rvec(xcopy
[ii
],state
->x
[ii
]);
2004 copy_rvec(vcopy
[ii
],state
->v
[ii
]);
2006 copy_mat(boxcopy
,state
->box
);
2011 /* for rerun MD always do Neighbour Searching */
2012 bNS
= (bFirstStep
|| ir
->nstlist
!= 0);
2017 /* Determine whether or not to do Neighbour Searching and LR */
2018 bNStList
= (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0);
2020 bNS
= (bFirstStep
|| bExchanged
|| bNStList
||
2021 (ir
->nstlist
== -1 && nlh
.nabnsb
> 0));
2023 if (bNS
&& ir
->nstlist
== -1)
2025 set_nlistheuristics(&nlh
,bFirstStep
|| bExchanged
,step
);
2029 if (terminate_now
> 0 || (terminate_now
< 0 && bNS
))
2034 /* Determine whether or not to update the Born radii if doing GB */
2035 bBornRadii
=bFirstStep
;
2036 if (ir
->implicit_solvent
&& (step
% ir
->nstgbradii
==0))
2041 do_log
= do_per_step(step
,ir
->nstlog
) || bFirstStep
|| bLastStep
;
2042 do_verbose
= bVerbose
&&
2043 (step
% stepout
== 0 || bFirstStep
|| bLastStep
);
2045 if (bNS
&& !(bFirstStep
&& ir
->bContinuation
&& !bRerunMD
))
2049 bMasterState
= TRUE
;
2053 bMasterState
= FALSE
;
2054 /* Correct the new box if it is too skewed */
2055 if (DYNAMIC_BOX(*ir
))
2057 if (correct_box(fplog
,step
,state
->box
,graph
))
2059 bMasterState
= TRUE
;
2062 if (DOMAINDECOMP(cr
) && bMasterState
)
2064 dd_collect_state(cr
->dd
,state
,state_global
);
2068 if (DOMAINDECOMP(cr
))
2070 /* Repartition the domain decomposition */
2071 wallcycle_start(wcycle
,ewcDOMDEC
);
2072 dd_partition_system(fplog
,step
,cr
,
2073 bMasterState
,nstglobalcomm
,
2074 state_global
,top_global
,ir
,
2075 state
,&f
,mdatoms
,top
,fr
,
2076 vsite
,shellfc
,constr
,
2077 nrnb
,wcycle
,do_verbose
);
2078 wallcycle_stop(wcycle
,ewcDOMDEC
);
2082 if (MASTER(cr
) && do_log
&& !bFFscan
)
2084 print_ebin_header(fplog
,step
,t
,state
->lambda
);
2087 if (ir
->efep
!= efepNO
)
2089 update_mdatoms(mdatoms
,state
->lambda
);
2092 if (bRerunMD
&& rerun_fr
.bV
)
2095 /* We need the kinetic energy at minus the half step for determining
2096 * the full step kinetic energy and possibly for T-coupling.*/
2097 /* This may not be quite working correctly yet . . . . */
2098 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2099 wcycle
,enerd
,NULL
,NULL
,NULL
,NULL
,mu_tot
,
2100 constr
,&chkpt
,&terminate
,&terminate_now
,&(nlh
.nabnsb
),state
->box
,
2101 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2102 (CGLO_RERUNMD
| CGLO_GSTAT
| CGLO_TEMPERATURE
) |
2103 (bEkinAveVel
?CGLO_EKINAVEVEL
:0));
2105 clear_mat(force_vir
);
2107 /* Ionize the atoms if necessary */
2110 ionize(fplog
,oenv
,mdatoms
,top_global
,t
,ir
,state
->x
,state
->v
,
2111 mdatoms
->start
,mdatoms
->start
+mdatoms
->homenr
,state
->box
,cr
);
2114 /* Update force field in ffscan program */
2117 if (update_forcefield(fplog
,
2119 mdatoms
->nr
,state
->x
,state
->box
)) {
2120 if (gmx_parallel_env())
2128 GMX_MPE_LOG(ev_timestep2
);
2130 if ((bNS
|| bLastStep
) && (step
> ir
->init_step
) && !bRerunMD
)
2132 bCPT
= (chkpt
> 0 || (bLastStep
&& (Flags
& MD_CONFOUT
)));
2143 /* Determine the pressure:
2144 * always when we want exact averages in the energy file,
2145 * at ns steps when we have pressure coupling,
2146 * otherwise only at energy output steps (set below).
2148 bNstEner
= (bGStatEveryStep
|| do_per_step(step
,ir
->nstcalcenergy
));
2149 bCalcEner
= bNstEner
;
2150 bCalcPres
= (bGStatEveryStep
|| (ir
->epc
!= epcNO
&& bNS
));
2152 /* Do we need global communication ? */
2153 bGStat
= (bCalcEner
|| bStopCM
||
2154 (ir
->nstlist
== -1 && !bRerunMD
&& step
>= nlh
.step_nscheck
));
2156 do_ene
= (do_per_step(step
,ir
->nstenergy
) || bLastStep
);
2158 if (do_ene
|| do_log
)
2165 /* these CGLO_ options remain the same throughout the iteration */
2166 cglo_flags
= ((bRerunMD
? CGLO_RERUNMD
: 0) |
2167 (bEkinAveVel
? CGLO_EKINAVEVEL
: 0) |
2168 (bStopCM
? CGLO_STOPCM
: 0) |
2169 (bGStat
? CGLO_GSTAT
: 0) |
2170 (bNEMD
? CGLO_NEMD
: 0));
2172 force_flags
= (GMX_FORCE_STATECHANGED
|
2173 ((DYNAMIC_BOX(*ir
) || bRerunMD
) ? GMX_FORCE_DYNAMICBOX
: 0) |
2174 GMX_FORCE_ALLFORCES
|
2175 (bNStList
? GMX_FORCE_DOLR
: 0) |
2177 (bCalcEner
? GMX_FORCE_VIRIAL
: 0) |
2178 (bDoDHDL
? GMX_FORCE_DHDL
: 0));
2182 /* Now is the time to relax the shells */
2183 count
=relax_shell_flexcon(fplog
,cr
,bVerbose
,bFFscan
? step
+1 : step
,
2185 bStopCM
,top
,top_global
,
2187 state
,f
,force_vir
,mdatoms
,
2188 nrnb
,wcycle
,graph
,groups
,
2189 shellfc
,fr
,bBornRadii
,t
,mu_tot
,
2190 state
->natoms
,&bConverged
,vsite
,
2201 /* The coordinates (x) are shifted (to get whole molecules)
2203 * This is parallellized as well, and does communication too.
2204 * Check comments in sim_util.c
2207 do_force(fplog
,cr
,ir
,step
,nrnb
,wcycle
,top
,top_global
,groups
,
2208 state
->box
,state
->x
,&state
->hist
,
2209 f
,force_vir
,mdatoms
,enerd
,fcd
,
2210 state
->lambda
,graph
,
2211 fr
,vsite
,mu_tot
,t
,fp_field
,ed
,bBornRadii
,
2212 (bNS
? GMX_FORCE_NS
: 0) | force_flags
);
2215 GMX_BARRIER(cr
->mpi_comm_mygroup
);
2219 mu_aver
= calc_mu_aver(cr
,state
->x
,mdatoms
->chargeA
,
2220 mu_tot
,&top_global
->mols
,mdatoms
,gnx
,grpindex
);
2223 if (bTCR
&& bFirstStep
)
2225 tcr
=init_coupling(fplog
,nfile
,fnm
,cr
,fr
,mdatoms
,&(top
->idef
));
2226 fprintf(fplog
,"Done init_coupling\n");
2230 /* ############### START FIRST UPDATE HALF-STEP ############### */
2232 if (!bStartingFromCpt
&& !bRerunMD
) {
2233 if (bVVAveVel
&& bInitStep
) {
2234 /* if using velocity verlet with full time step Ekin, take the first half step only to compute the
2235 virial for the first step. From there, revert back to the initial coordinates
2236 so that the input is actually the initial step */
2237 copy_rvecn(state
->v
,cbuf
,0,state
->natoms
); /* should make this better for parallelizing? */
2240 /* this is for NHC in the Ekin(t+dt/2) version of vv */
2241 if (!bInitStep
|| !bEkinAveVel
)
2243 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[1],bEkinAveVel
);
2246 update_coords(fplog
,step
,ir
,mdatoms
,state
,
2247 f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2248 ekind
,M
,wcycle
,upd
,bInitStep
,etrtVELOCITY
,cr
,nrnb
,constr
,&top
->idef
);
2250 bIterateFirstHalf
= bIterate
&& (!bInitStep
);
2251 bFirstIterate
= TRUE
;
2252 /* for iterations, we save these vectors, as we will be self-consistently iterating
2254 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
2256 /* save the state */
2257 if (bIterateFirstHalf
) {
2258 copy_coupling_state(state
,bufstate
,ekind
,ekind_save
);
2261 while (bIterateFirstHalf
|| bFirstIterate
)
2263 if (bIterateFirstHalf
)
2265 copy_coupling_state(bufstate
,state
,ekind_save
,ekind
);
2267 if (bIterateFirstHalf
)
2269 if (bFirstIterate
&& bTrotter
)
2271 /* The first time, we need a decent first estimate
2272 of veta(t+dt) to compute the constraints. Do
2273 this by computing the box volume part of the
2274 trotter integration at this time. Nothing else
2275 should be changed by this routine here. If
2276 !(first time), we start with the previous value
2279 veta_save
= state
->veta
;
2280 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[0],FALSE
);
2281 vetanew
= state
->veta
;
2282 state
->veta
= veta_save
;
2287 if ( !bRerunMD
|| rerun_fr
.bV
|| bForceUpdate
) { /* Why is rerun_fr.bV here? Unclear. */
2288 wallcycle_start(wcycle
,ewcUPDATE
);
2291 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
2292 &top
->idef
,shake_vir
,NULL
,
2293 cr
,nrnb
,wcycle
,upd
,constr
,
2294 bInitStep
,TRUE
,bCalcPres
,vetanew
);
2296 if (!bOK
&& !bFFscan
)
2298 gmx_fatal(FARGS
,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2303 { /* Need to unshift here if a do_force has been
2304 called in the previous step */
2305 unshift_self(graph
,state
->box
,state
->x
);
2308 /* if bEkinAveVel, then compute ekin and shake_vir, otherwise, just compute shake_vir */
2309 /* first step, we need to compute the virial */
2311 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2312 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
2313 constr
,&chkpt
,&terminate
,&terminate_now
,&(nlh
.nabnsb
),state
->box
,
2314 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2315 (cglo_flags
| CGLO_ENERGY
) |
2316 ((bIterateFirstHalf
| bInitStep
| IR_NPT_TROTTER(ir
)) ? (cglo_flags
| CGLO_PRESSURE
):0) |
2317 (((bEkinAveVel
&&(!bInitStep
)) || (!bEkinAveVel
&& IR_NPT_TROTTER(ir
)))? (cglo_flags
| CGLO_TEMPERATURE
):0) |
2318 ((bEkinAveVel
|| (!bEkinAveVel
&& IR_NPT_TROTTER(ir
))) ? CGLO_EKINAVEVEL
: 0) |
2319 (bIterate
? CGLO_ITERATE
: 0) |
2320 (bFirstIterate
? CGLO_FIRSTITERATE
: 0) |
2321 (cglo_flags
| CGLO_SCALEEKIN
)
2324 /* explanation of above:
2325 a) We compute Ekin at the full time step
2326 if 1) we are using the AveVel Ekin, and it's not the
2327 initial step, or 2) if we are using AveEkin, but need the full
2328 time step kinetic energy for the pressure.
2329 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
2330 EkinAveVel because it's needed for the pressure */
2332 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
2335 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[2],bEkinAveVel
);
2338 if (done_iterating(cr
,fplog
,&bFirstIterate
,&bIterateFirstHalf
,state
->veta
,&vetanew
,0)) break;
2341 if (bTrotter
&& !bInitStep
) {
2342 copy_mat(shake_vir
,state
->vir_prev
);
2343 if (IR_NVT_TROTTER(ir
) && bEkinAveVel
) {
2344 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
2345 enerd
->term
[F_TEMP
] = sum_ekin(&(ir
->opts
),ekind
,NULL
,bEkinAveVel
,FALSE
,FALSE
);
2346 enerd
->term
[F_EKIN
] = trace(ekind
->ekin
);
2349 /* if it's the initial step, we performed this first step just to get the constraint virial */
2350 if (bInitStep
&& bVVAveVel
) {
2351 copy_rvecn(cbuf
,state
->v
,0,state
->natoms
);
2354 if (fr
->bSepDVDL
&& fplog
&& do_log
)
2356 fprintf(fplog
,sepdvdlformat
,"Constraint",0.0,dvdl
);
2358 enerd
->term
[F_DHDL_CON
] += dvdl
;
2360 wallcycle_stop(wcycle
,ewcUPDATE
);
2361 GMX_MPE_LOG(ev_timestep1
);
2363 /* MRS -- done iterating -- compute the conserved quantity */
2367 if (IR_NVT_TROTTER(ir
) || IR_NPT_TROTTER(ir
))
2370 NPT_energy(ir
,state
->nosehoover_xi
,state
->nosehoover_vxi
,state
->veta
, state
->box
, &MassQ
);
2371 if ((ir
->eDispCorr
!= edispcEnerPres
) && (ir
->eDispCorr
!= edispcAllEnerPres
))
2373 last_conserved
-= enerd
->term
[F_DISPCORR
];
2377 last_ekin
= enerd
->term
[F_EKIN
]; /* does this get preserved through checkpointing? */
2381 /* ######## END FIRST UPDATE STEP ############## */
2382 /* ######## If doing VV, we now have v(dt) ###### */
2384 /* ################## START TRAJECTORY OUTPUT ################# */
2386 /* Now we have the energies and forces corresponding to the
2387 * coordinates at time t. We must output all of this before
2389 * for RerunMD t is read from input trajectory
2391 GMX_MPE_LOG(ev_output_start
);
2393 bX
= do_per_step(step
,ir
->nstxout
);
2394 bV
= do_per_step(step
,ir
->nstvout
);
2395 bF
= do_per_step(step
,ir
->nstfout
);
2396 bXTC
= do_per_step(step
,ir
->nstxtcout
);
2400 fcReportProgress( ir
->nsteps
, step
);
2402 bX
= bX
|| bLastStep
; /*enforce writing positions and velocities
2404 bV
= bV
|| bLastStep
;
2406 int nthreads
=(cr
->nthreads
==0 ? 1 : cr
->nthreads
);
2407 int nnodes
=(cr
->nnodes
==0 ? 1 : cr
->nnodes
);
2410 /*Gromacs drives checkpointing; no ||
2411 fcCheckPointPendingThreads(cr->nodeid,
2413 /* sync bCPT and fc record-keeping */
2414 if (bCPT
&& MASTER(cr
))
2415 fcRequestCheckPoint();
2420 if (bX
|| bV
|| bF
|| bXTC
|| bCPT
)
2422 wallcycle_start(wcycle
,ewcTRAJ
);
2425 if (state
->flags
& (1<<estLD_RNG
))
2427 get_stochd_state(upd
,state
);
2433 state_global
->ekinstate
.bUpToDate
= FALSE
;
2437 update_ekinstate(&state_global
->ekinstate
,ekind
);
2438 state_global
->ekinstate
.bUpToDate
= TRUE
;
2440 update_energyhistory(&state_global
->enerhist
,mdebin
);
2443 write_traj(fplog
,cr
,fp_trn
,bX
,bV
,bF
,fp_xtc
,bXTC
,ir
->xtcprec
,
2444 fn_cpt
,bCPT
,top_global
,ir
->eI
,ir
->simulation_part
,
2445 step
,t
,state
,state_global
,f
,f_global
,&n_xtc
,&x_xtc
);
2447 if (bLastStep
&& step_rel
== ir
->nsteps
&&
2448 (Flags
& MD_CONFOUT
) && MASTER(cr
) &&
2449 !bRerunMD
&& !bFFscan
)
2451 /* x and v have been collected in write_traj */
2452 fprintf(stderr
,"\nWriting final coordinates.\n");
2453 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
&&
2456 /* Make molecules whole only for confout writing */
2457 do_pbc_mtop(fplog
,ir
->ePBC
,state
->box
,top_global
,state_global
->x
);
2459 write_sto_conf_mtop(ftp2fn(efSTO
,nfile
,fnm
),
2460 *top_global
->name
,top_global
,
2461 state_global
->x
,state_global
->v
,
2462 ir
->ePBC
,state
->box
);
2465 wallcycle_stop(wcycle
,ewcTRAJ
);
2467 GMX_MPE_LOG(ev_output_finish
);
2469 /* kludge -- virial is lost with restart for NPT control. Must restart */
2470 if (bStartingFromCpt
&& bVV
)
2472 copy_mat(state
->vir_prev
,shake_vir
);
2474 /* ################## END TRAJECTORY OUTPUT ################ */
2476 /* Determine the pressure:
2477 * always when we want exact averages in the energy file,
2478 * at ns steps when we have pressure coupling,
2479 * otherwise only at energy output steps (set below).
2482 bNstEner
= (bGStatEveryStep
|| do_per_step(step
,ir
->nstcalcenergy
));
2483 bCalcEner
= bNstEner
;
2484 bCalcPres
= (bGStatEveryStep
|| (ir
->epc
!= epcNO
&& bNS
));
2486 /* Do we need global communication ? */
2487 bGStat
= (bGStatEveryStep
|| bStopCM
|| bNS
||
2488 (ir
->nstlist
== -1 && !bRerunMD
&& step
>= nlh
.step_nscheck
));
2490 do_ene
= (do_per_step(step
,ir
->nstenergy
) || bLastStep
);
2492 if (do_ene
|| do_log
)
2498 bFirstIterate
= TRUE
;
2500 /* for iterations, we save these vectors, as we will be redoing the calculations */
2503 copy_coupling_state(state
,bufstate
,ekind
,ekind_save
);
2506 while (bIterate
|| bFirstIterate
)
2508 /* We now restore these vectors to redo the calculation with improved extended variables */
2511 copy_coupling_state(bufstate
,state
,ekind_save
,ekind
);
2514 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
2515 so scroll down for that logic */
2517 /* ######### START SECOND UPDATE STEP ################# */
2519 GMX_MPE_LOG(ev_update_start
);
2521 if (!bRerunMD
|| rerun_fr
.bV
|| bForceUpdate
)
2523 wallcycle_start(wcycle
,ewcUPDATE
);
2527 /* Box is changed in update() when we do pressure coupling,
2528 * but we should still use the old box for energy corrections and when
2529 * writing it to the energy file, so it matches the trajectory files for
2530 * the same timestep above. Make a copy in a separate array.
2533 copy_mat(state
->box
,lastbox
);
2535 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS*/
2540 tracevir
= trace(shake_vir
);
2544 /* we use a new value of scalevir to converge the iterations faster */
2545 scalevir
= tracevir
/trace(shake_vir
);
2546 msmul(shake_vir
,scalevir
,shake_vir
);
2547 m_add(force_vir
,shake_vir
,total_vir
);
2548 clear_mat(shake_vir
);
2550 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[3],bEkinAveVel
);
2552 /* We can only do Berendsen coupling after we have summed the kinetic
2553 * energy or virial. Since the happens in global_state after update,
2554 * we should only do it at step % nstlist = 1 with bGStatEveryStep=FALSE.
2557 update_extended(fplog
,step
,ir
,mdatoms
,state
,ekind
,pcoupl_mu
,M
,wcycle
,
2558 upd
,bInitStep
,FALSE
,&MassQ
);
2560 /* velocity (for VV) */
2561 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2562 ekind
,M
,wcycle
,upd
,FALSE
,etrtVELOCITY
,cr
,nrnb
,constr
,&top
->idef
);
2564 /* above, initialize just copies ekinh into ekin, it doesn't copy
2565 /* position (for VV), and entire integrator for MD */
2569 copy_rvecn(state
->x
,cbuf
,0,state
->natoms
);
2572 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2573 ekind
,M
,wcycle
,upd
,bInitStep
,etrtPOSITION
,cr
,nrnb
,constr
,&top
->idef
);
2575 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
2576 &top
->idef
,shake_vir
,force_vir
,
2577 cr
,nrnb
,wcycle
,upd
,constr
,
2578 bInitStep
,FALSE
,bCalcPres
,state
->veta
);
2582 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2583 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
2584 constr
,&chkpt
,&terminate
,&terminate_now
,&(nlh
.nabnsb
),lastbox
,
2585 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2586 (cglo_flags
| CGLO_TEMPERATURE
) |
2587 (cglo_flags
& ~CGLO_PRESSURE
) |
2588 (cglo_flags
& ~CGLO_ENERGY
)
2590 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[4],bEkinAveVel
);
2591 copy_rvecn(cbuf
,state
->x
,0,state
->natoms
);
2592 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2593 ekind
,M
,wcycle
,upd
,bInitStep
,etrtPOSITION
,cr
,nrnb
,constr
,&top
->idef
);
2594 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
2595 /* are the small terms in the shake_vir here due
2596 * to numerical errors, or are they important
2597 * physically? I'm thinking they are just errors, but not completely sure. */
2598 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
2599 &top
->idef
,tmp_vir
,force_vir
,
2600 cr
,nrnb
,wcycle
,upd
,constr
,
2601 bInitStep
,FALSE
,bCalcPres
,state
->veta
);
2602 m_add(shake_vir
,tmp_vir
,shake_vir
);
2605 if (!bOK
&& !bFFscan
)
2607 gmx_fatal(FARGS
,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2610 if (fr
->bSepDVDL
&& fplog
&& do_log
)
2612 fprintf(fplog
,sepdvdlformat
,"Constraint",0.0,dvdl
);
2614 enerd
->term
[F_DHDL_CON
] += dvdl
;
2615 wallcycle_stop(wcycle
,ewcUPDATE
);
2619 /* Need to unshift here */
2620 unshift_self(graph
,state
->box
,state
->x
);
2623 GMX_BARRIER(cr
->mpi_comm_mygroup
);
2624 GMX_MPE_LOG(ev_update_finish
);
2628 wallcycle_start(wcycle
,ewcVSITECONSTR
);
2631 shift_self(graph
,state
->box
,state
->x
);
2633 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,state
->v
,
2634 top
->idef
.iparams
,top
->idef
.il
,
2635 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
2639 unshift_self(graph
,state
->box
,state
->x
);
2641 wallcycle_stop(wcycle
,ewcVSITECONSTR
);
2644 /* ############## IF NOT VV, CALCULATE EKIN HERE - ALWAYS CALCULATE PRESSURE ############ */
2645 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2646 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
2647 constr
,&chkpt
,&terminate
,&terminate_now
,&(nlh
.nabnsb
),lastbox
,
2648 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2649 (!bVV
? (cglo_flags
| CGLO_ENERGY
):0) |
2650 (cglo_flags
| CGLO_PRESSURE
) |
2651 (!bVV
? (cglo_flags
| CGLO_TEMPERATURE
):0) |
2652 (bIterate
? CGLO_ITERATE
: 0) |
2653 (bFirstIterate
? CGLO_FIRSTITERATE
: 0)
2655 /* bIterate is set to keep it from eliminating the old ekinh */
2656 /* ############# END CALC EKIN AND PRESSURE ################# */
2658 if (done_iterating(cr
,fplog
,&bFirstIterate
,&bIterate
,trace(shake_vir
),&tracevir
,1)) break;
2661 update_box(fplog
,step
,ir
,mdatoms
,state
,graph
,f
,
2662 ir
->nstlist
==-1 ? &nlh
.scale_tot
: NULL
,pcoupl_mu
,nrnb
,wcycle
,upd
,bInitStep
,FALSE
);
2664 /* ################# END UPDATE STEP 2 ################# */
2665 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
2667 /* Determine the wallclock run time up till now */
2668 run_time
= (double)time(NULL
) - (double)runtime
->real
;
2670 /* Check whether everything is still allright */
2671 if ((bGotTermSignal
|| bGotUsr1Signal
) && !bHandledSignal
)
2673 if (bGotTermSignal
|| ir
->nstlist
== 0)
2683 terminate_now
= terminate
;
2688 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2689 bGotTermSignal
? "TERM" : "USR1",
2690 terminate
==-1 ? "NS " : "");
2694 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2695 bGotTermSignal
? "TERM" : "USR1",
2696 terminate
==-1 ? "NS " : "");
2698 bHandledSignal
=TRUE
;
2700 else if (MASTER(cr
) && (bNS
|| ir
->nstlist
<= 0) &&
2701 (max_hours
> 0 && run_time
> max_hours
*60.0*60.0*0.99) &&
2704 /* Signal to terminate the run */
2705 terminate
= (ir
->nstlist
== 0 ? 1 : -1);
2708 terminate_now
= terminate
;
2712 fprintf(fplog
,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step
,sbuf
),max_hours
*0.99);
2714 fprintf(stderr
, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step
,sbuf
),max_hours
*0.99);
2717 if (ir
->nstlist
== -1 && !bRerunMD
)
2719 /* When bGStatEveryStep=FALSE, global_stat is only called
2720 * when we check the atom displacements, not at NS steps.
2721 * This means that also the bonded interaction count check is not
2722 * performed immediately after NS. Therefore a few MD steps could
2723 * be performed with missing interactions.
2724 * But wrong energies are never written to file,
2725 * since energies are only written after global_stat
2728 if (step
>= nlh
.step_nscheck
)
2730 nlh
.nabnsb
= natoms_beyond_ns_buffer(ir
,fr
,&top
->cgs
,
2731 nlh
.scale_tot
,state
->x
);
2735 /* This is not necessarily true,
2736 * but step_nscheck is determined quite conservatively.
2742 /* In parallel we only have to check for checkpointing in steps
2743 * where we do global communication,
2744 * otherwise the other nodes don't know.
2746 if (MASTER(cr
) && ((bGStat
|| !PAR(cr
)) &&
2749 run_time
>= nchkpt
*cpt_period
*60.0)))
2758 /* The coordinates (x) were unshifted in update */
2759 if (bFFscan
&& (shellfc
==NULL
|| bConverged
))
2761 if (print_forcefield(fplog
,enerd
->term
,mdatoms
->homenr
,
2763 &(top_global
->mols
),mdatoms
->massT
,pres
))
2765 if (gmx_parallel_env()) {
2768 fprintf(stderr
,"\n");
2774 /* We will not sum ekinh_old,
2775 * so signal that we still have to do it.
2777 bSumEkinhOld
= TRUE
;
2782 /* Only do GCT when the relaxation of shells (minimization) has converged,
2783 * otherwise we might be coupling to bogus energies.
2784 * In parallel we must always do this, because the other sims might
2788 /* Since this is called with the new coordinates state->x, I assume
2789 * we want the new box state->box too. / EL 20040121
2791 do_coupling(fplog
,oenv
,nfile
,fnm
,tcr
,t
,step
,enerd
->term
,fr
,
2793 mdatoms
,&(top
->idef
),mu_aver
,
2794 top_global
->mols
.nr
,cr
,
2795 state
->box
,total_vir
,pres
,
2796 mu_tot
,state
->x
,f
,bConverged
);
2800 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
2802 sum_dhdl(enerd
,state
->lambda
,ir
);
2803 /* use the directly determined last velocity, not actually the averaged half steps */
2804 if (bTrotter
&& bEkinAveVel
) {
2805 enerd
->term
[F_EKIN
] = last_ekin
;
2807 enerd
->term
[F_ETOT
] = enerd
->term
[F_EPOT
] + enerd
->term
[F_EKIN
];
2816 if (IR_NVT_TROTTER(ir
)) {
2817 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] + last_conserved
;
2819 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] +
2820 NPT_energy(ir
,state
->nosehoover_xi
,
2821 state
->nosehoover_vxi
,state
->veta
,state
->box
,&MassQ
);
2825 enerd
->term
[F_ECONSERVED
] =
2826 enerd
->term
[F_ETOT
] + vrescale_energy(&(ir
->opts
),
2827 state
->therm_integral
);
2833 /* Check for excessively large energies */
2837 real etot_max
= 1e200
;
2839 real etot_max
= 1e30
;
2841 if (fabs(enerd
->term
[F_ETOT
]) > etot_max
)
2843 fprintf(stderr
,"Energy too large (%g), giving up\n",
2844 enerd
->term
[F_ETOT
]);
2847 /* ######### END PREPARING EDR OUTPUT ########### */
2849 /* Time for performance */
2850 if (((step
% stepout
) == 0) || bLastStep
)
2852 runtime_upd_proc(runtime
);
2862 upd_mdebin(mdebin
,bDoDHDL
? fp_dhdl
: NULL
,TRUE
,
2863 t
,mdatoms
->tmass
,enerd
,state
,lastbox
,
2864 shake_vir
,force_vir
,total_vir
,pres
,
2865 ekind
,mu_tot
,constr
);
2869 upd_mdebin_step(mdebin
);
2872 do_dr
= do_per_step(step
,ir
->nstdisreout
);
2873 do_or
= do_per_step(step
,ir
->nstorireout
);
2875 print_ebin(fp_ene
,do_ene
,do_dr
,do_or
,do_log
?fplog
:NULL
,step
,t
,
2876 eprNORMAL
,bCompact
,mdebin
,fcd
,groups
,&(ir
->opts
));
2878 if (ir
->ePull
!= epullNO
)
2880 pull_print_output(ir
->pull
,step
,t
);
2883 if (do_per_step(step
,ir
->nstlog
))
2885 if(fflush(fplog
) != 0)
2887 gmx_fatal(FARGS
,"Cannot flush logfile - maybe you are out of quota?");
2893 /* Remaining runtime */
2894 if (MULTIMASTER(cr
) && do_verbose
)
2898 fprintf(stderr
,"\n");
2900 print_time(stderr
,runtime
,step
,ir
);
2903 /* Replica exchange */
2905 if ((repl_ex_nst
> 0) && (step
> 0) && !bLastStep
&&
2906 do_per_step(step
,repl_ex_nst
))
2908 bExchanged
= replica_exchange(fplog
,cr
,repl_ex
,
2909 state_global
,enerd
->term
,
2912 if (bExchanged
&& PAR(cr
))
2914 if (DOMAINDECOMP(cr
))
2916 dd_partition_system(fplog
,step
,cr
,TRUE
,1,
2917 state_global
,top_global
,ir
,
2918 state
,&f
,mdatoms
,top
,fr
,
2919 vsite
,shellfc
,constr
,
2924 bcast_state(cr
,state
,FALSE
);
2930 bStartingFromCpt
= FALSE
;
2932 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
2933 /* Complicated conditional when bGStatEveryStep=FALSE.
2934 * We can not just use bGStat, since then the simulation results
2935 * would depend on nstenergy and nstlog or step_nscheck.
2937 if (((state
->flags
& (1<<estPRES_PREV
)) || (state
->flags
& (1<<estVIR_PREV
))) &&
2939 (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0) ||
2940 (ir
->nstlist
< 0 && nlh
.nabnsb
> 0) ||
2941 (ir
->nstlist
== 0 && bGStat
)))
2943 /* Store the pressure in t_state for pressure coupling
2944 * at the next MD step.
2946 if (state
->flags
& (1<<estPRES_PREV
))
2948 copy_mat(pres
,state
->pres_prev
);
2952 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
2956 /* read next frame from input trajectory */
2957 bNotLastFrame
= read_next_frame(oenv
,status
,&rerun_fr
);
2960 if (!bRerunMD
|| !rerun_fr
.bStep
)
2962 /* increase the MD step number */
2967 cycles
= wallcycle_stop(wcycle
,ewcSTEP
);
2968 if (DOMAINDECOMP(cr
) && wcycle
)
2970 dd_cycles_add(cr
->dd
,cycles
,ddCyclStep
);
2973 if (step_rel
== wcycle_get_reset_counters(wcycle
))
2975 /* Reset all the counters related to performance over the run */
2976 reset_all_counters(fplog
,cr
,step
,&step_rel
,ir
,wcycle
,nrnb
,runtime
);
2977 wcycle_set_reset_counters(wcycle
, 0);
2980 /* End of main MD loop */
2984 runtime_end(runtime
);
2991 if (!(cr
->duty
& DUTY_PME
))
2993 /* Tell the PME only node to finish */
2999 if (bGStatEveryStep
&& !bRerunMD
)
3001 print_ebin(fp_ene
,FALSE
,FALSE
,FALSE
,fplog
,step
,t
,
3002 eprAVER
,FALSE
,mdebin
,fcd
,groups
,&(ir
->opts
));
3012 gmx_fio_fclose(fp_dhdl
);
3016 gmx_fio_fclose(fp_field
);
3021 if (ir
->nstlist
== -1 && nlh
.nns
> 0 && fplog
)
3023 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
)));
3024 fprintf(fplog
,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh
.ab
/nlh
.nns
);
3027 if (shellfc
&& fplog
)
3029 fprintf(fplog
,"Fraction of iterations that converged: %.2f %%\n",
3030 (nconverged
*100.0)/step_rel
);
3031 fprintf(fplog
,"Average number of force evaluations per MD step: %.2f\n\n",
3035 if (repl_ex_nst
> 0 && MASTER(cr
))
3037 print_replica_exchange_statistics(fplog
,repl_ex
);
3040 runtime
->nsteps_done
= step_rel
;