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 real ekin
,temp
,prescorr
,enercorr
,dvdlcorr
;
193 /* translate CGLO flags to booleans */
194 bRerunMD
= flags
& CGLO_RERUNMD
;
195 bEkinAveVel
= flags
& CGLO_EKINAVEVEL
;
196 bStopCM
= flags
& CGLO_STOPCM
;
197 bGStat
= flags
& CGLO_GSTAT
;
198 bNEMD
= flags
& CGLO_NEMD
;
199 bIterate
= flags
& CGLO_ITERATE
;
200 bFirstIterate
= flags
& CGLO_FIRSTITERATE
;
201 bCopyEkinh
= flags
& CGLO_COPYEKINH
;
202 bEner
= flags
& CGLO_ENERGY
;
203 bTemp
= flags
& CGLO_TEMPERATURE
;
204 bPres
= flags
& CGLO_PRESSURE
;
205 bReadEkin
= flags
& CGLO_READEKIN
;
206 bScaleEkin
= flags
& CGLO_SCALEEKIN
;
208 /* decide when to calculate temperature and pressure. */
210 /* in initalization, it sums the shake virial in vv, and to
211 sums ekinh_old in leapfrog (or if we are calculating ekinh_old for other reasons */
213 bVV
= (ir
->eI
== eiVV
);
215 /* ########## Kinetic energy ############## */
219 /* Non-equilibrium MD: this is parallellized, but only does communication
220 * when there really is NEMD.
223 if (PAR(cr
) && (bNEMD
))
225 accumulate_u(cr
,&(ir
->opts
),ekind
);
230 restore_ekinstate_from_state(cr
,ekind
,&state_global
->ekinstate
);
234 calc_ke_part(state
,&(ir
->opts
),mdatoms
,ekind
,nrnb
,bEkinAveVel
,bIterate
);
239 /* Calculate center of mass velocity if necessary, also parallellized */
240 if (bStopCM
&& !bRerunMD
&& bEner
)
242 calc_vcm_grp(fplog
,mdatoms
->start
,mdatoms
->homenr
,mdatoms
,
243 state
->x
,state
->v
,vcm
);
247 if (bTemp
|| bPres
|| bEner
)
251 /* We will not sum ekinh_old,
252 * so signal that we still have to do it.
254 *bSumEkinhOld
= TRUE
;
260 GMX_MPE_LOG(ev_global_stat_start
);
261 global_stat(fplog
,gstat
,cr
,enerd
,force_vir
,shake_vir
,mu_tot
,
262 ir
,ekind
,constr
,vcm
,NULL
,NULL
,terminate
,top_global
,state
,
263 *bSumEkinhOld
,flags
);
264 GMX_MPE_LOG(ev_global_stat_finish
);
269 terminate_now
= terminate
;
272 wallcycle_stop(wcycle
,ewcMoveE
);
273 *bSumEkinhOld
= FALSE
;
277 if (!bNEMD
&& debug
&& bTemp
&& (vcm
->nr
> 0))
280 mdatoms
->start
,mdatoms
->start
+mdatoms
->homenr
,
281 state
->v
,vcm
->group_p
[0],
282 mdatoms
->massT
,mdatoms
->tmass
,ekind
->ekin
);
286 /* Do center of mass motion removal */
287 if (bStopCM
&& !bRerunMD
&& !bFirstHalf
) /* fix this! */
289 check_cm_grp(fplog
,vcm
,ir
,1);
290 do_stopcm_grp(fplog
,mdatoms
->start
,mdatoms
->homenr
,mdatoms
->cVCM
,
291 state
->x
,state
->v
,vcm
);
292 inc_nrnb(nrnb
,eNR_STOPCM
,mdatoms
->homenr
);
298 /* Sum the kinetic energies of the groups & calc temp */
299 enerd
->term
[F_TEMP
] = sum_ekin(&(ir
->opts
),ekind
,&(enerd
->term
[F_DKDL
]),
300 bEkinAveVel
||bReadEkin
,bIterate
,bScaleEkin
);
301 /* three main: VV with AveVel, vv with AveEkin, leap with AveEkin. Leap with AveVel is also
302 an option, but not supported now. Additionally, if we are doing iterations.
303 bCopyHalf: if TRUE, we simply copy the ekinh directly to ekin, multiplying be ekinscale_nhc.
304 bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale.
305 If FALSE, we average ekinh_old and ekinh*ekinscale_nhc.
306 bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
307 If FALSE, we go ahead and erase.
309 enerd
->term
[F_EKIN
] = trace(ekind
->ekin
);
312 /* ########## Long range energy information ###### */
316 calc_dispcorr(fplog
,ir
,fr
,0,top_global
,box
,state
->lambda
,
317 corr_pres
,corr_vir
,&prescorr
,&enercorr
,&dvdlcorr
);
320 if (bEner
&& bFirstIterate
)
322 enerd
->term
[F_DISPCORR
] = enercorr
;
323 enerd
->term
[F_EPOT
] += enercorr
;
324 enerd
->term
[F_DVDL
] += dvdlcorr
;
325 if (fr
->efep
!= efepNO
) {
326 enerd
->dvdl_lin
+= dvdlcorr
;
330 /* ########## Now pressure ############## */
334 m_add(force_vir
,shake_vir
,total_vir
);
336 /* Calculate pressure and apply LR correction if PPPM is used.
337 * Use the box from last timestep since we already called update().
340 enerd
->term
[F_PRES
] = calc_pres(fr
->ePBC
,ir
->nwall
,box
,ekind
->ekin
,total_vir
,pres
,
341 (fr
->eeltype
==eelPPPM
)?enerd
->term
[F_COUL_RECIP
]:0.0);
343 /* Calculate long range corrections to pressure and energy */
344 /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
345 and computes enerd->term[F_DISPCORR]. Also modifies the
346 total_vir and pres tesors */
348 m_add(total_vir
,corr_vir
,total_vir
);
349 m_add(pres
,corr_pres
,pres
);
350 enerd
->term
[F_PDISPCORR
] = prescorr
;
351 enerd
->term
[F_PRES
] += prescorr
;
352 *pcurr
= enerd
->term
[F_PRES
];
353 /* calculate temperature using virial */
354 enerd
->term
[F_VTEMP
] = calc_temp(trace(total_vir
),ir
->opts
.nrdf
[0]);
360 struct mdrunner_arglist
374 const char *dddlb_opt
;
387 int ret
; /* return value */
391 static void mdrunner_start_fn(void *arg
)
393 struct mdrunner_arglist
*mda
=(struct mdrunner_arglist
*)arg
;
394 struct mdrunner_arglist mc
=*mda
; /* copy the arg list to make sure
395 that it's thread-local. This doesn't
396 copy pointed-to items, of course,
397 but those are all const. */
398 t_commrec
*cr
; /* we need a local version of this */
400 t_filenm
*fnm
=dup_tfn(mc
.nfile
, mc
.fnm
);
402 cr
=init_par_threads(mc
.cr
);
409 mda
->ret
=mdrunner(fplog
, cr
, mc
.nfile
, mc
.fnm
, mc
.oenv
, mc
.bVerbose
,
410 mc
.bCompact
, mc
.nstglobalcomm
,
411 mc
.ddxyz
, mc
.dd_node_order
, mc
.rdd
,
412 mc
.rconstr
, mc
.dddlb_opt
, mc
.dlb_scale
,
413 mc
.ddcsx
, mc
.ddcsy
, mc
.ddcsz
, mc
.nstepout
, mc
.nmultisim
,
414 mc
.repl_ex_nst
, mc
.repl_ex_seed
, mc
.pforce
,
415 mc
.cpt_period
, mc
.max_hours
, mc
.Flags
);
420 int mdrunner_threads(int nthreads
,
421 FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
422 const output_env_t oenv
, bool bVerbose
,bool bCompact
,
424 ivec ddxyz
,int dd_node_order
,real rdd
,real rconstr
,
425 const char *dddlb_opt
,real dlb_scale
,
426 const char *ddcsx
,const char *ddcsy
,const char *ddcsz
,
427 int nstepout
,int nmultisim
,int repl_ex_nst
,
428 int repl_ex_seed
, real pforce
,real cpt_period
,
429 real max_hours
, unsigned long Flags
)
432 /* first check whether we even need to start tMPI */
435 ret
=mdrunner(fplog
, cr
, nfile
, fnm
, oenv
, bVerbose
, bCompact
,
437 ddxyz
, dd_node_order
, rdd
, rconstr
, dddlb_opt
, dlb_scale
,
438 ddcsx
, ddcsy
, ddcsz
, nstepout
, nmultisim
, repl_ex_nst
,
439 repl_ex_seed
, pforce
, cpt_period
, max_hours
, Flags
);
444 struct mdrunner_arglist mda
;
445 /* fill the data structure to pass as void pointer to thread start fn */
451 mda
.bVerbose
=bVerbose
;
452 mda
.bCompact
=bCompact
;
453 mda
.nstglobalcomm
=nstglobalcomm
;
454 mda
.ddxyz
[XX
]=ddxyz
[XX
];
455 mda
.ddxyz
[YY
]=ddxyz
[YY
];
456 mda
.ddxyz
[ZZ
]=ddxyz
[ZZ
];
457 mda
.dd_node_order
=dd_node_order
;
460 mda
.dddlb_opt
=dddlb_opt
;
461 mda
.dlb_scale
=dlb_scale
;
465 mda
.nstepout
=nstepout
;
466 mda
.nmultisim
=nmultisim
;
467 mda
.repl_ex_nst
=repl_ex_nst
;
468 mda
.repl_ex_seed
=repl_ex_seed
;
470 mda
.cpt_period
=cpt_period
;
471 mda
.max_hours
=max_hours
;
474 fprintf(stderr
, "Starting %d threads\n",nthreads
);
476 tMPI_Init_fn(nthreads
, mdrunner_start_fn
, (void*)(&mda
) );
480 gmx_comm("Multiple threads requested but not compiled with threads");
487 int mdrunner(FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
488 const output_env_t oenv
, bool bVerbose
,bool bCompact
,
490 ivec ddxyz
,int dd_node_order
,real rdd
,real rconstr
,
491 const char *dddlb_opt
,real dlb_scale
,
492 const char *ddcsx
,const char *ddcsy
,const char *ddcsz
,
493 int nstepout
,int nmultisim
,int repl_ex_nst
,int repl_ex_seed
,
494 real pforce
,real cpt_period
,real max_hours
,
497 double nodetime
=0,realtime
;
498 t_inputrec
*inputrec
;
505 gmx_mtop_t
*mtop
=NULL
;
506 t_mdatoms
*mdatoms
=NULL
;
510 gmx_pme_t
*pmedata
=NULL
;
511 gmx_vsite_t
*vsite
=NULL
;
513 int i
,m
,nChargePerturbed
=-1,status
,nalloc
;
515 gmx_wallcycle_t wcycle
;
516 bool bReadRNG
,bReadEkin
;
518 gmx_runtime_t runtime
;
520 gmx_large_int_t reset_counters
;
523 /* Essential dynamics */
524 if (opt2bSet("-ei",nfile
,fnm
))
526 /* Open input and output files, allocate space for ED data structure */
527 ed
= ed_open(nfile
,fnm
,cr
);
535 if (bVerbose
&& SIMMASTER(cr
))
537 fprintf(stderr
,"Getting Loaded...\n");
540 if (Flags
& MD_APPENDFILES
)
547 /* The master thread on the master node reads from disk,
548 * then distributes everything to the other processors.
551 list
= (SIMMASTER(cr
) && !(Flags
& MD_APPENDFILES
)) ? (LIST_SCALARS
| LIST_INPUTREC
) : 0;
554 init_parallel(fplog
, opt2fn_master("-s",nfile
,fnm
,cr
),cr
,
555 inputrec
,mtop
,state
,list
);
560 /* Read a file for a single processor */
562 init_single(fplog
,inputrec
,ftp2fn(efTPX
,nfile
,fnm
),mtop
,state
);
564 if (!EEL_PME(inputrec
->coulombtype
) || (Flags
& MD_PARTDEC
))
570 fcRegisterSteps(inputrec
->nsteps
,inputrec
->init_step
);
573 /* NMR restraints must be initialized before load_checkpoint,
574 * since with time averaging the history is added to t_state.
575 * For proper consistency check we therefore need to extend
577 * So the PME-only nodes (if present) will also initialize
578 * the distance restraints.
582 /* This needs to be called before read_checkpoint to extend the state */
583 init_disres(fplog
,mtop
,inputrec
,cr
,Flags
& MD_PARTDEC
,fcd
,state
);
585 if (gmx_mtop_ftype_count(mtop
,F_ORIRES
) > 0)
587 if (PAR(cr
) && !(Flags
& MD_PARTDEC
))
589 gmx_fatal(FARGS
,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
591 /* Orientation restraints */
594 init_orires(fplog
,mtop
,state
->x
,inputrec
,cr
->ms
,&(fcd
->orires
),
599 if (DEFORM(*inputrec
))
601 /* Store the deform reference box before reading the checkpoint */
604 copy_mat(state
->box
,box
);
608 gmx_bcast(sizeof(box
),box
,cr
);
610 /* Because we do not have the update struct available yet
611 * in which the reference values should be stored,
612 * we store them temporarily in static variables.
613 * This should be thread safe, since they are only written once
614 * and with identical values.
617 tMPI_Thread_mutex_lock(&box_mutex
);
619 init_step_tpx
= inputrec
->init_step
;
620 copy_mat(box
,box_tpx
);
622 tMPI_Thread_mutex_unlock(&box_mutex
);
626 if (opt2bSet("-cpi",nfile
,fnm
))
628 /* Check if checkpoint file exists before doing continuation.
629 * This way we can use identical input options for the first and subsequent runs...
631 if( gmx_fexist_master(opt2fn_master("-cpi",nfile
,fnm
,cr
),cr
) )
633 load_checkpoint(opt2fn_master("-cpi",nfile
,fnm
,cr
),fplog
,
634 cr
,Flags
& MD_PARTDEC
,ddxyz
,
635 inputrec
,state
,&bReadRNG
,&bReadEkin
,
636 (Flags
& MD_APPENDFILES
));
640 Flags
|= MD_READ_RNG
;
644 Flags
|= MD_READ_EKIN
;
649 if (MASTER(cr
) && (Flags
& MD_APPENDFILES
))
651 fplog
= gmx_log_open(ftp2fn(efLOG
,nfile
,fnm
),cr
,!(Flags
& MD_SEPPOT
),
657 copy_mat(state
->box
,box
);
662 gmx_bcast(sizeof(box
),box
,cr
);
665 if (bVerbose
&& SIMMASTER(cr
))
667 fprintf(stderr
,"Loaded with Money\n\n");
670 if (PAR(cr
) && !((Flags
& MD_PARTDEC
) || EI_TPI(inputrec
->eI
)))
672 cr
->dd
= init_domain_decomposition(fplog
,cr
,Flags
,ddxyz
,rdd
,rconstr
,
679 make_dd_communicators(fplog
,cr
,dd_node_order
);
681 /* Set overallocation to avoid frequent reallocation of arrays */
682 set_over_alloc_dd(TRUE
);
686 cr
->duty
= (DUTY_PP
| DUTY_PME
);
687 npme_major
= cr
->nnodes
;
689 if (inputrec
->ePBC
== epbcSCREW
)
692 "pbc=%s is only implemented with domain decomposition",
693 epbc_names
[inputrec
->ePBC
]);
699 /* After possible communicator splitting in make_dd_communicators.
700 * we can set up the intra/inter node communication.
702 gmx_setup_nodecomm(fplog
,cr
);
705 wcycle
= wallcycle_init(fplog
,cr
);
708 /* Master synchronizes its value of reset_counters with all nodes
709 * including PME only nodes */
710 reset_counters
= wcycle_get_reset_counters(wcycle
);
711 gmx_bcast_sim(sizeof(reset_counters
),&reset_counters
,cr
);
712 wcycle_set_reset_counters(wcycle
, reset_counters
);
717 if (cr
->duty
& DUTY_PP
)
719 /* For domain decomposition we allocate dynamically
720 * in dd_partition_system.
722 if (DOMAINDECOMP(cr
))
724 bcast_state_setup(cr
,state
);
734 bcast_state(cr
,state
,TRUE
);
738 /* Dihedral Restraints */
739 if (gmx_mtop_ftype_count(mtop
,F_DIHRES
) > 0)
741 init_dihres(fplog
,mtop
,inputrec
,fcd
);
744 /* Initiate forcerecord */
746 init_forcerec(fplog
,oenv
,fr
,fcd
,inputrec
,mtop
,cr
,box
,FALSE
,
747 opt2fn("-table",nfile
,fnm
),
748 opt2fn("-tablep",nfile
,fnm
),
749 opt2fn("-tableb",nfile
,fnm
),FALSE
,pforce
);
751 /* version for PCA_NOT_READ_NODE (see md.c) */
752 /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
753 "nofile","nofile","nofile",FALSE,pforce);
755 fr
->bSepDVDL
= ((Flags
& MD_SEPPOT
) == MD_SEPPOT
);
757 /* Initialize QM-MM */
760 init_QMMMrec(cr
,box
,mtop
,inputrec
,fr
);
763 /* Initialize the mdatoms structure.
764 * mdatoms is not filled with atom data,
765 * as this can not be done now with domain decomposition.
767 mdatoms
= init_mdatoms(fplog
,mtop
,inputrec
->efep
!=efepNO
);
769 /* Initialize the virtual site communication */
770 vsite
= init_vsite(mtop
,cr
);
772 calc_shifts(box
,fr
->shift_vec
);
774 /* With periodic molecules the charge groups should be whole at start up
775 * and the virtual sites should not be far from their proper positions.
777 if (!inputrec
->bContinuation
&& MASTER(cr
) &&
778 !(inputrec
->ePBC
!= epbcNONE
&& inputrec
->bPeriodicMols
))
780 /* Make molecules whole at start of run */
781 if (fr
->ePBC
!= epbcNONE
)
783 do_pbc_first_mtop(fplog
,inputrec
->ePBC
,box
,mtop
,state
->x
);
787 /* Correct initial vsite positions are required
788 * for the initial distribution in the domain decomposition
789 * and for the initial shell prediction.
791 construct_vsites_mtop(fplog
,vsite
,mtop
,state
->x
);
795 /* Initiate PPPM if necessary */
796 if (fr
->eeltype
== eelPPPM
)
798 if (mdatoms
->nChargePerturbed
)
800 gmx_fatal(FARGS
,"Free energy with %s is not implemented",
801 eel_names
[fr
->eeltype
]);
803 status
= gmx_pppm_init(fplog
,cr
,oenv
,FALSE
,TRUE
,box
,
804 getenv("GMXGHAT"),inputrec
, (Flags
& MD_REPRODUCIBLE
));
807 gmx_fatal(FARGS
,"Error %d initializing PPPM",status
);
811 if (EEL_PME(fr
->eeltype
))
813 ewaldcoeff
= fr
->ewaldcoeff
;
814 pmedata
= &fr
->pmedata
;
823 /* This is a PME only node */
825 /* We don't need the state */
828 ewaldcoeff
= calc_ewaldcoeff(inputrec
->rcoulomb
, inputrec
->ewald_rtol
);
832 /* Initiate PME if necessary,
833 * either on all nodes or on dedicated PME nodes only. */
834 if (EEL_PME(inputrec
->coulombtype
))
838 nChargePerturbed
= mdatoms
->nChargePerturbed
;
840 if (cr
->npmenodes
> 0)
842 /* The PME only nodes need to know nChargePerturbed */
843 gmx_bcast_sim(sizeof(nChargePerturbed
),&nChargePerturbed
,cr
);
845 if (cr
->duty
& DUTY_PME
)
847 status
= gmx_pme_init(pmedata
,cr
,npme_major
,inputrec
,
848 mtop
? mtop
->natoms
: 0,nChargePerturbed
,
849 (Flags
& MD_REPRODUCIBLE
));
852 gmx_fatal(FARGS
,"Error %d initializing PME",status
);
858 if (integrator
[inputrec
->eI
].func
== do_md
)
860 /* Turn on signal handling on all nodes */
862 * (A user signal from the PME nodes (if any)
863 * is communicated to the PP nodes.
865 if (getenv("GMX_NO_TERM") == NULL
)
869 fprintf(debug
,"Installing signal handler for SIGTERM\n");
871 signal(SIGTERM
,signal_handler
);
874 if (getenv("GMX_NO_USR1") == NULL
)
878 fprintf(debug
,"Installing signal handler for SIGUSR1\n");
880 signal(SIGUSR1
,signal_handler
);
885 if (cr
->duty
& DUTY_PP
)
887 if (inputrec
->ePull
!= epullNO
)
889 /* Initialize pull code */
890 init_pull(fplog
,inputrec
,nfile
,fnm
,mtop
,cr
,oenv
,
891 EI_DYNAMICS(inputrec
->eI
) && MASTER(cr
),Flags
);
894 constr
= init_constraints(fplog
,mtop
,inputrec
,ed
,state
,cr
);
896 if (DOMAINDECOMP(cr
))
898 dd_init_bondeds(fplog
,cr
->dd
,mtop
,vsite
,constr
,inputrec
,
899 Flags
& MD_DDBONDCHECK
,fr
->cginfo_mb
);
901 set_dd_parameters(fplog
,cr
->dd
,dlb_scale
,inputrec
,fr
,&ddbox
);
903 setup_dd_grid(fplog
,cr
->dd
);
906 /* Now do whatever the user wants us to do (how flexible...) */
907 integrator
[inputrec
->eI
].func(fplog
,cr
,nfile
,fnm
,
908 oenv
,bVerbose
,bCompact
,
911 nstepout
,inputrec
,mtop
,
913 mdatoms
,nrnb
,wcycle
,ed
,fr
,
914 repl_ex_nst
,repl_ex_seed
,
915 cpt_period
,max_hours
,
919 if (inputrec
->ePull
!= epullNO
)
921 finish_pull(fplog
,inputrec
->pull
);
927 gmx_pmeonly(*pmedata
,cr
,nrnb
,wcycle
,ewaldcoeff
,FALSE
,inputrec
);
930 if (EI_DYNAMICS(inputrec
->eI
) || EI_TPI(inputrec
->eI
))
932 /* Some timing stats */
935 if (runtime
.proc
== 0)
937 runtime
.proc
= runtime
.real
;
946 wallcycle_stop(wcycle
,ewcRUN
);
948 /* Finish up, write some stuff
949 * if rerunMD, don't write last frame again
951 finish_run(fplog
,cr
,ftp2fn(efSTO
,nfile
,fnm
),
952 inputrec
,nrnb
,wcycle
,&runtime
,
953 EI_DYNAMICS(inputrec
->eI
) && !MULTISIM(cr
));
955 /* Does what it says */
956 print_date_and_time(fplog
,cr
->nodeid
,"Finished mdrun",&runtime
);
958 /* Close logfile already here if we were appending to it */
959 if (MASTER(cr
) && (Flags
& MD_APPENDFILES
))
961 gmx_log_close(fplog
);
968 else if(bGotUsr1Signal
)
980 static void md_print_warning(const t_commrec
*cr
,FILE *fplog
,const char *buf
)
984 fprintf(stderr
,"\n%s\n",buf
);
988 fprintf(fplog
,"\n%s\n",buf
);
992 static bool done_iterating(const t_commrec
*cr
,FILE *fplog
, bool *bFirstIterate
, bool *bIterate
, real fom
, real
*newf
, int n
)
997 /* monitor convergence, and use a secant search to propose new
1000 The secant method computes x_{i+1} = x_{i} - f(x_{i}) * ---------------------
1001 f(x_{i}) - f(x_{i-1})
1003 The function we are trying to zero is fom-x, where fom is the
1004 "figure of merit" which is the pressure (or the veta value) we
1005 would get by putting in an old value of the pressure or veta into
1006 the incrementor function for the step or half step. I have
1007 verified that this gives the same answer as self consistent
1008 iteration, usually in many fewer steps, especially for small tau_p.
1010 We could possibly eliminate an iteration with proper use
1011 of the value from the previous step, but that would take a bit
1012 more bookkeeping, especially for veta, since tests indicate the
1013 function of veta on the last step is not sufficiently close to
1014 guarantee convergence this step. This is
1015 good enough for now. On my tests, I could use tau_p down to
1016 0.02, which is smaller that would ever be necessary in
1017 practice. Generally, 3-5 iterations will be sufficient */
1019 /* Definitions for convergence. Could be optimized . . . */
1021 #define CONVERGEITER 0.000000001
1023 #define CONVERGEITER 0.0001
1025 #define MAXITERCONST 200
1027 /* used to escape out of cyclic traps because of limited numberical precision */
1029 #define FALLBACK 1000
1031 static real f
,fprev
,x
,xprev
;
1033 static real allrelerr
[MAXITERCONST
+2];
1041 *bFirstIterate
= FALSE
;
1049 f
= fom
-x
; /* we want to zero this difference */
1050 if ((iter_i
> 1) && (iter_i
< MAXITERCONST
))
1058 *newf
= x
- (x
-xprev
)*(f
)/(f
-fprev
);
1063 /* just use self-consistent iteration the first step to initialize, or
1064 if it's not converging (which happens occasionally -- need to investigate why) */
1068 /* Consider a slight shortcut allowing us to exit one sooner -- we check the
1069 difference between the closest of x and xprev to the new
1070 value. To be 100% certain, we should check the difference between
1071 the last result, and the previous result, or
1073 relerr = (fabs((x-xprev)/fom));
1075 but this is pretty much never necessary under typical conditions.
1076 Checking numerically, it seems to lead to almost exactly the same
1077 trajectories, but there are small differences out a few decimal
1078 places in the pressure, and eventually in the v_eta, but it could
1081 if (fabs(*newf-x) < fabs(*newf - xprev)) { xmin = x;} else { xmin = xprev;}
1082 relerr = (fabs((*newf-xmin) / *newf));
1085 relerr
= (fabs((f
-fprev
)/fom
));
1087 allrelerr
[iter_i
] = relerr
;
1093 fprintf(debug
,"Iterating NPT constraints #%i: %6i %20.12f%14.6g%20.12f\n",n
,iter_i
,fom
,relerr
,*newf
);
1096 if ((relerr
< CONVERGEITER
) || (fom
==0))
1101 fprintf(debug
,"Iterating NPT constraints #%i: CONVERGED\n",n
);
1105 if (iter_i
> MAXITERCONST
)
1108 for (i
=0;i
<CYCLEMAX
;i
++) {
1109 if (allrelerr
[(MAXITERCONST
-2*CYCLEMAX
)-i
] == allrelerr
[iter_i
-1]) {
1111 if (relerr
> CONVERGEITER
*FALLBACK
) {incycle
= FALSE
; break;}
1116 /* we are trapped in a numerical attractor, and can't converge any more, and are close to the final result.
1117 Better to give up here and proceed with a warning than have the simulation die.
1119 sprintf(buf
,"numerical convergence issues with NPT, relative error only %10.5g, continuing\n",relerr
);
1120 md_print_warning(cr
,fplog
,buf
);
1123 gmx_fatal(FARGS
,"Could not converge NPT constraints\n");
1136 static void check_nst_param(FILE *fplog
,t_commrec
*cr
,
1137 const char *desc_nst
,int nst
,
1138 const char *desc_p
,int *p
)
1142 if (*p
> 0 && *p
% nst
!= 0)
1144 /* Round up to the next multiple of nst */
1145 *p
= ((*p
)/nst
+ 1)*nst
;
1146 sprintf(buf
,"NOTE: %s changes %s to %d\n",desc_nst
,desc_p
,*p
);
1147 md_print_warning(cr
,fplog
,buf
);
1151 static void reset_all_counters(FILE *fplog
,t_commrec
*cr
,
1152 gmx_large_int_t step
,
1153 gmx_large_int_t
*step_rel
,t_inputrec
*ir
,
1154 gmx_wallcycle_t wcycle
,t_nrnb
*nrnb
,
1155 gmx_runtime_t
*runtime
)
1157 char buf
[STRLEN
],sbuf
[22];
1159 /* Reset all the counters related to performance over the run */
1160 sprintf(buf
,"Step %s: resetting all time and cycle counters\n",
1161 gmx_step_str(step
,sbuf
));
1162 md_print_warning(cr
,fplog
,buf
);
1164 wallcycle_stop(wcycle
,ewcRUN
);
1165 wallcycle_reset_all(wcycle
);
1166 if (DOMAINDECOMP(cr
))
1168 reset_dd_statistics_counters(cr
->dd
);
1171 ir
->init_step
+= *step_rel
;
1172 ir
->nsteps
-= *step_rel
;
1174 wallcycle_start(wcycle
,ewcRUN
);
1175 runtime_start(runtime
);
1176 print_date_and_time(fplog
,cr
->nodeid
,"Restarted time",runtime
);
1179 static int check_nstglobalcomm(FILE *fplog
,t_commrec
*cr
,
1180 int nstglobalcomm
,t_inputrec
*ir
)
1184 if (!EI_DYNAMICS(ir
->eI
))
1189 if (nstglobalcomm
== -1)
1191 if (ir
->nstcalcenergy
== 0 && ir
->nstlist
== 0)
1194 if (ir
->nstenergy
> 0 && ir
->nstenergy
< nstglobalcomm
)
1196 nstglobalcomm
= ir
->nstenergy
;
1201 /* We assume that if nstcalcenergy > nstlist,
1202 * nstcalcenergy is a multiple of nstlist.
1204 if (ir
->nstcalcenergy
== 0 ||
1205 (ir
->nstlist
> 0 && ir
->nstlist
< ir
->nstcalcenergy
))
1207 nstglobalcomm
= ir
->nstlist
;
1211 nstglobalcomm
= ir
->nstcalcenergy
;
1217 if (ir
->nstlist
> 0 &&
1218 nstglobalcomm
> ir
->nstlist
&& nstglobalcomm
% ir
->nstlist
!= 0)
1220 nstglobalcomm
= (nstglobalcomm
/ ir
->nstlist
)*ir
->nstlist
;
1221 sprintf(buf
,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm
);
1222 md_print_warning(cr
,fplog
,buf
);
1224 if (nstglobalcomm
> ir
->nstcalcenergy
)
1226 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
1227 "nstcalcenergy",&ir
->nstcalcenergy
);
1230 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
1231 "nstenergy",&ir
->nstenergy
);
1233 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
1234 "nstlog",&ir
->nstlog
);
1237 if (ir
->comm_mode
!= ecmNO
&& ir
->nstcomm
< nstglobalcomm
)
1239 sprintf(buf
,"WARNING: Changing nstcomm from %d to %d\n",
1240 ir
->nstcomm
,nstglobalcomm
);
1241 md_print_warning(cr
,fplog
,buf
);
1242 ir
->nstcomm
= nstglobalcomm
;
1245 return nstglobalcomm
;
1248 static void check_ir_old_tpx_versions(t_commrec
*cr
,FILE *fplog
,
1249 t_inputrec
*ir
,gmx_mtop_t
*mtop
)
1251 /* Check required for old tpx files */
1252 if (IR_TWINRANGE(*ir
) && ir
->nstlist
> 1 &&
1253 ir
->nstcalcenergy
% ir
->nstlist
!= 0)
1255 md_print_warning(cr
,fplog
,"Old tpr file with twin-range settings: modifiying energy calculation and/or T/P-coupling frequencies");
1257 if (gmx_mtop_ftype_count(mtop
,F_CONSTR
) +
1258 gmx_mtop_ftype_count(mtop
,F_CONSTRNC
) > 0 &&
1259 ir
->eConstrAlg
== econtSHAKE
)
1261 md_print_warning(cr
,fplog
,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect");
1262 if (ir
->epc
!= epcNO
)
1264 gmx_fatal(FARGS
,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
1267 check_nst_param(fplog
,cr
,"nstlist",ir
->nstlist
,
1268 "nstcalcenergy",&ir
->nstcalcenergy
);
1270 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
1271 "nstenergy",&ir
->nstenergy
);
1272 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
1273 "nstlog",&ir
->nstlog
);
1274 if (ir
->efep
!= efepNO
)
1276 check_nst_param(fplog
,cr
,"nstdhdl",ir
->nstdhdl
,
1277 "nstenergy",&ir
->nstenergy
);
1282 bool bGStatEveryStep
;
1283 gmx_large_int_t step_ns
;
1284 gmx_large_int_t step_nscheck
;
1285 gmx_large_int_t nns
;
1295 static void reset_nlistheuristics(gmx_nlheur_t
*nlh
,gmx_large_int_t step
)
1299 nlh
->step_nscheck
= step
;
1302 static void init_nlistheuristics(gmx_nlheur_t
*nlh
,
1303 bool bGStatEveryStep
,gmx_large_int_t step
)
1305 nlh
->bGStatEveryStep
= bGStatEveryStep
;
1312 reset_nlistheuristics(nlh
,step
);
1315 static void update_nliststatistics(gmx_nlheur_t
*nlh
,gmx_large_int_t step
)
1317 gmx_large_int_t nl_lt
;
1318 char sbuf
[22],sbuf2
[22];
1320 /* Determine the neighbor list life time */
1321 nl_lt
= step
- nlh
->step_ns
;
1324 fprintf(debug
,"%d atoms beyond ns buffer, updating neighbor list after %s steps\n",nlh
->nabnsb
,gmx_step_str(nl_lt
,sbuf
));
1328 nlh
->s2
+= nl_lt
*nl_lt
;
1329 nlh
->ab
+= nlh
->nabnsb
;
1330 if (nlh
->lt_runav
== 0)
1332 nlh
->lt_runav
= nl_lt
;
1333 /* Initialize the fluctuation average
1334 * such that at startup we check after 0 steps.
1336 nlh
->lt_runav2
= sqr(nl_lt
/2.0);
1338 /* Running average with 0.9 gives an exp. history of 9.5 */
1339 nlh
->lt_runav2
= 0.9*nlh
->lt_runav2
+ 0.1*sqr(nlh
->lt_runav
- nl_lt
);
1340 nlh
->lt_runav
= 0.9*nlh
->lt_runav
+ 0.1*nl_lt
;
1341 if (nlh
->bGStatEveryStep
)
1343 /* Always check the nlist validity */
1344 nlh
->step_nscheck
= step
;
1348 /* We check after: <life time> - 2*sigma
1349 * The factor 2 is quite conservative,
1350 * but we assume that with nstlist=-1 the user
1351 * prefers exact integration over performance.
1353 nlh
->step_nscheck
= step
1354 + (int)(nlh
->lt_runav
- 2.0*sqrt(nlh
->lt_runav2
)) - 1;
1358 fprintf(debug
,"nlist life time %s run av. %4.1f sig %3.1f check %s check with -gcom %d\n",
1359 gmx_step_str(nl_lt
,sbuf
),nlh
->lt_runav
,sqrt(nlh
->lt_runav2
),
1360 gmx_step_str(nlh
->step_nscheck
-step
+1,sbuf2
),
1361 (int)(nlh
->lt_runav
- 2.0*sqrt(nlh
->lt_runav2
)));
1365 static void set_nlistheuristics(gmx_nlheur_t
*nlh
,bool bReset
,gmx_large_int_t step
)
1371 reset_nlistheuristics(nlh
,step
);
1375 update_nliststatistics(nlh
,step
);
1378 nlh
->step_ns
= step
;
1379 /* Initialize the cumulative coordinate scaling matrix */
1380 clear_mat(nlh
->scale_tot
);
1381 for(d
=0; d
<DIM
; d
++)
1383 nlh
->scale_tot
[d
][d
] = 1.0;
1387 double do_md(FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
1388 const output_env_t oenv
, bool bVerbose
,bool bCompact
,
1390 gmx_vsite_t
*vsite
,gmx_constr_t constr
,
1391 int stepout
,t_inputrec
*ir
,
1392 gmx_mtop_t
*top_global
,
1394 t_state
*state_global
,
1396 t_nrnb
*nrnb
,gmx_wallcycle_t wcycle
,
1397 gmx_edsam_t ed
,t_forcerec
*fr
,
1398 int repl_ex_nst
,int repl_ex_seed
,
1399 real cpt_period
,real max_hours
,
1400 unsigned long Flags
,
1401 gmx_runtime_t
*runtime
)
1403 int fp_trn
=0,fp_xtc
=0;
1404 ener_file_t fp_ene
=NULL
;
1405 gmx_large_int_t step
,step_rel
;
1407 FILE *fp_dhdl
=NULL
,*fp_field
=NULL
;
1410 bool bGStatEveryStep
,bGStat
,bNstEner
,bCalcPres
,bCalcEner
;
1411 bool bNS
,bNStList
,bSimAnn
,bStopCM
,bRerunMD
,bNotLastFrame
=FALSE
,
1412 bFirstStep
,bStateFromTPX
,bInitStep
,bLastStep
,bEkinAveVel
,
1413 bBornRadii
,bStartingFromCpt
,bVVAveVel
,bVVAveEkin
;
1415 bool bNEMD
,do_ene
,do_log
,do_verbose
,bRerunWarnNoV
=TRUE
,
1416 bForceUpdate
=FALSE
,bX
,bV
,bF
,bXTC
,bCPT
;
1418 int force_flags
,cglo_flags
;
1419 tensor force_vir
,shake_vir
,total_vir
,tmp_vir
,pres
;
1424 matrix
*scale_tot
,pcoupl_mu
,M
,ebox
;
1426 t_trxframe rerun_fr
;
1427 gmx_repl_ex_t repl_ex
=NULL
;
1429 /* Booleans (disguised as a reals) to checkpoint and terminate mdrun */
1430 real chkpt
=0,terminate
=0,terminate_now
=0;
1432 gmx_localtop_t
*top
;
1433 t_mdebin
*mdebin
=NULL
;
1434 t_state
*state
=NULL
;
1435 rvec
*f_global
=NULL
;
1438 gmx_enerdata_t
*enerd
;
1440 gmx_global_stat_t gstat
;
1441 gmx_update_t upd
=NULL
;
1442 t_graph
*graph
=NULL
;
1445 gmx_groups_t
*groups
;
1446 gmx_ekindata_t
*ekind
, *ekind_save
;
1447 gmx_shellfc_t shellfc
;
1448 int count
,nconverged
=0;
1452 bool bTCR
=FALSE
,bConverged
=TRUE
,bOK
,bSumEkinhOld
,bExchanged
;
1454 bool bVV
,bIterate
,bIterateFirstHalf
,bFirstIterate
,bTemp
,bPres
,bTrotter
;
1455 real temp0
,mu_aver
=0,dvdl
;
1457 atom_id
*grpindex
=NULL
;
1459 t_coupl_rec
*tcr
=NULL
;
1460 rvec
*xcopy
=NULL
,*vcopy
=NULL
,*cbuf
=NULL
;
1461 matrix boxcopy
,lastbox
;
1463 real fom
,oldfom
,veta_save
,pcurr
,scalevir
,tracevir
;
1466 int reset_counters
=-1;
1467 real last_conserved
= 0;
1472 char sbuf
[22],sbuf2
[22];
1473 bool bHandledSignal
=FALSE
;
1475 /* Temporary addition for FAHCORE checkpointing */
1479 /* Check for special mdrun options */
1480 bRerunMD
= (Flags
& MD_RERUN
);
1481 bIonize
= (Flags
& MD_IONIZE
);
1482 bFFscan
= (Flags
& MD_FFSCAN
);
1483 bAppend
= (Flags
& MD_APPENDFILES
);
1484 /* md-vv uses averaged half step velocities to determine temperature unless defined otherwise by GMX_EKIN_AVE_EKIN;
1485 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
1486 bVV
= (ir
->eI
==eiVV
);
1488 bEkinAveVel
= (getenv("GMX_EKIN_AVE_EKIN")==NULL
);
1490 bEkinAveVel
= (getenv("GMX_EKIN_AVE_VEL")!=NULL
);
1492 bVVAveVel
= bVV
&& bEkinAveVel
;
1493 bVVAveEkin
= bVV
&& !bEkinAveVel
;
1495 if (bVV
) /* to store the initial velocities while computing virial */
1497 snew(cbuf
,top_global
->natoms
);
1499 /* all the iteratative cases - only if there are constraints */
1500 bIterate
= ((IR_NPT_TROTTER(ir
)) && (constr
) && (!bRerunMD
));
1501 bTrotter
= (bVV
&& (IR_NPT_TROTTER(ir
) || (IR_NVT_TROTTER(ir
))));
1505 /* Since we don't know if the frames read are related in any way,
1506 * rebuild the neighborlist at every step.
1509 ir
->nstcalcenergy
= 1;
1513 check_ir_old_tpx_versions(cr
,fplog
,ir
,top_global
);
1515 nstglobalcomm
= check_nstglobalcomm(fplog
,cr
,nstglobalcomm
,ir
);
1516 bGStatEveryStep
= (nstglobalcomm
== 1);
1518 if (!bGStatEveryStep
&& ir
->nstlist
== -1 && fplog
!= NULL
)
1521 "To reduce the energy communication with nstlist = -1\n"
1522 "the neighbor list validity should not be checked at every step,\n"
1523 "this means that exact integration is not guaranteed.\n"
1524 "The neighbor list validity is checked after:\n"
1525 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
1526 "In most cases this will result in exact integration.\n"
1527 "This reduces the energy communication by a factor of 2 to 3.\n"
1528 "If you want less energy communication, set nstlist > 3.\n\n");
1531 if (bRerunMD
|| bFFscan
)
1535 groups
= &top_global
->groups
;
1537 /* Initial values */
1538 init_md(fplog
,cr
,ir
,oenv
,&t
,&t0
,&state_global
->lambda
,&lam0
,
1539 nrnb
,top_global
,&upd
,
1540 nfile
,fnm
,&fp_trn
,&fp_xtc
,&fp_ene
,&fn_cpt
,
1541 &fp_dhdl
,&fp_field
,&mdebin
,
1542 force_vir
,shake_vir
,mu_tot
,&bNEMD
,&bSimAnn
,&vcm
,Flags
);
1544 clear_mat(total_vir
);
1546 /* Energy terms and groups */
1548 init_enerdata(top_global
->groups
.grps
[egcENER
].nr
,ir
->n_flambda
,enerd
);
1549 if (DOMAINDECOMP(cr
))
1555 snew(f
,top_global
->natoms
);
1558 /* Kinetic energy data */
1560 init_ekindata(fplog
,top_global
,&(ir
->opts
),ekind
);
1561 /* needed for iteration of constraints */
1563 init_ekindata(fplog
,top_global
,&(ir
->opts
),ekind_save
);
1564 /* Copy the cos acceleration to the groups struct */
1565 ekind
->cosacc
.cos_accel
= ir
->cos_accel
;
1567 gstat
= global_stat_init(ir
);
1570 /* Check for polarizable models and flexible constraints */
1571 shellfc
= init_shell_flexcon(fplog
,
1572 top_global
,n_flexible_constraints(constr
),
1573 (ir
->bContinuation
||
1574 (DOMAINDECOMP(cr
) && !MASTER(cr
))) ?
1575 NULL
: state_global
->x
);
1580 tMPI_Thread_mutex_lock(&box_mutex
);
1582 set_deform_reference_box(upd
,init_step_tpx
,box_tpx
);
1584 tMPI_Thread_mutex_unlock(&box_mutex
);
1589 double io
= compute_io(ir
,top_global
->natoms
,groups
,mdebin
->ebin
->nener
,1);
1590 if ((io
> 2000) && MASTER(cr
))
1592 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
1596 if (DOMAINDECOMP(cr
)) {
1597 top
= dd_init_local_top(top_global
);
1600 dd_init_local_state(cr
->dd
,state_global
,state
);
1602 if (DDMASTER(cr
->dd
) && ir
->nstfout
) {
1603 snew(f_global
,state_global
->natoms
);
1607 /* Initialize the particle decomposition and split the topology */
1608 top
= split_system(fplog
,top_global
,ir
,cr
);
1610 pd_cg_range(cr
,&fr
->cg0
,&fr
->hcg
);
1611 pd_at_range(cr
,&a0
,&a1
);
1613 top
= gmx_mtop_generate_local_top(top_global
,ir
);
1616 a1
= top_global
->natoms
;
1619 state
= partdec_init_local_state(cr
,state_global
);
1622 atoms2md(top_global
,ir
,0,NULL
,a0
,a1
-a0
,mdatoms
);
1625 set_vsite_top(vsite
,top
,mdatoms
,cr
);
1628 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
) {
1629 graph
= mk_graph(fplog
,&(top
->idef
),0,top_global
->natoms
,FALSE
,FALSE
);
1633 make_local_shells(cr
,mdatoms
,shellfc
);
1636 if (ir
->pull
&& PAR(cr
)) {
1637 dd_make_local_pull_groups(NULL
,ir
->pull
,mdatoms
);
1641 if (DOMAINDECOMP(cr
))
1643 /* Distribute the charge groups over the nodes from the master node */
1644 dd_partition_system(fplog
,ir
->init_step
,cr
,TRUE
,1,
1645 state_global
,top_global
,ir
,
1646 state
,&f
,mdatoms
,top
,fr
,
1647 vsite
,shellfc
,constr
,
1651 /* If not DD, copy gb data */
1652 if(ir
->implicit_solvent
&& !DOMAINDECOMP(cr
))
1654 make_local_gb(cr
,fr
->born
,ir
->gb_algorithm
);
1657 update_mdatoms(mdatoms
,state
->lambda
);
1661 /* Update mdebin with energy history if appending to output files */
1662 if ( Flags
& MD_APPENDFILES
)
1664 restore_energyhistory_from_state(mdebin
,&state_global
->enerhist
);
1666 /* Set the initial energy history in state to zero by updating once */
1667 update_energyhistory(&state_global
->enerhist
,mdebin
);
1670 if ((state
->flags
& (1<<estLD_RNG
)) && (Flags
& MD_READ_RNG
)) {
1671 /* Set the random state if we read a checkpoint file */
1672 set_stochd_state(upd
,state
);
1675 /* Initialize constraints */
1677 if (!DOMAINDECOMP(cr
))
1678 set_constraints(constr
,top
,ir
,mdatoms
,cr
);
1681 /* Check whether we have to GCT stuff */
1682 bTCR
= ftp2bSet(efGCT
,nfile
,fnm
);
1685 fprintf(stderr
,"Will do General Coupling Theory!\n");
1687 gnx
= top_global
->mols
.nr
;
1689 for(i
=0; (i
<gnx
); i
++) {
1694 if (repl_ex_nst
> 0 && MASTER(cr
))
1695 repl_ex
= init_replica_exchange(fplog
,cr
->ms
,state_global
,ir
,
1696 repl_ex_nst
,repl_ex_seed
);
1698 if (!ir
->bContinuation
&& !bRerunMD
)
1700 if (mdatoms
->cFREEZE
&& (state
->flags
& (1<<estV
)))
1702 /* Set the velocities of frozen particles to zero */
1703 for(i
=mdatoms
->start
; i
<mdatoms
->start
+mdatoms
->homenr
; i
++)
1705 for(m
=0; m
<DIM
; m
++)
1707 if (ir
->opts
.nFreeze
[mdatoms
->cFREEZE
[i
]][m
])
1717 /* Constrain the initial coordinates and velocities */
1718 do_constrain_first(fplog
,constr
,ir
,mdatoms
,state
,f
,
1719 graph
,cr
,nrnb
,fr
,top
,shake_vir
);
1723 /* Construct the virtual sites for the initial configuration */
1724 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,NULL
,
1725 top
->idef
.iparams
,top
->idef
.il
,
1726 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
1732 /* would need to modify if vv starts with full time step */
1733 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
1734 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
1735 constr
,&chkpt
,&terminate
,&terminate_now
,&(nlh
.nabnsb
),state
->box
,
1736 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
1738 (bVV
? CGLO_PRESSURE
:0) |
1739 (CGLO_TEMPERATURE
) |
1740 (bRerunMD
? CGLO_RERUNMD
:0) |
1741 (bEkinAveVel
? CGLO_EKINAVEVEL
:0) |
1743 ((Flags
& MD_READ_EKIN
) ? CGLO_READEKIN
:0)));
1744 /* I'm assuming we need global communication the first time */
1745 /* Calculate the initial half step temperature, and save the ekinh_old */
1746 if (!(Flags
& MD_STARTFROMCPT
))
1748 for(i
=0; (i
<ir
->opts
.ngtc
); i
++)
1750 copy_mat(ekind
->tcstat
[i
].ekinh
,ekind
->tcstat
[i
].ekinh_old
);
1753 temp0
= enerd
->term
[F_TEMP
];
1755 /* Initiate data for the special cases */
1758 bufstate
= init_bufstate(state
->natoms
,(ir
->opts
.ngtc
+1)*NNHCHAIN
); /* extra state for barostat */
1763 snew(xcopy
,state
->natoms
);
1764 snew(vcopy
,state
->natoms
);
1765 copy_rvecn(state
->x
,xcopy
,0,state
->natoms
);
1766 copy_rvecn(state
->v
,vcopy
,0,state
->natoms
);
1767 copy_mat(state
->box
,boxcopy
);
1770 /* need to make an initial call to get the Trotter variables set. */
1771 trotter_seq
= init_trotter(ir
,state
,&MassQ
,bTrotter
);
1775 if (constr
&& !ir
->bContinuation
&& ir
->eConstrAlg
== econtLINCS
)
1778 "RMS relative constraint deviation after constraining: %.2e\n",
1779 constr_rmsd(constr
,FALSE
));
1781 fprintf(fplog
,"Initial temperature: %g K\n",enerd
->term
[F_TEMP
]);
1784 fprintf(stderr
,"starting md rerun '%s', reading coordinates from"
1785 " input trajectory '%s'\n\n",
1786 *(top_global
->name
),opt2fn("-rerun",nfile
,fnm
));
1789 fprintf(stderr
,"Calculated time to finish depends on nsteps from "
1790 "run input file,\nwhich may not correspond to the time "
1791 "needed to process input trajectory.\n\n");
1796 fprintf(stderr
,"starting mdrun '%s'\n",
1797 *(top_global
->name
));
1798 if (ir
->init_step
> 0)
1800 fprintf(stderr
,"%s steps, %8.1f ps (continuing from step %s, %8.1f ps).\n",
1801 gmx_step_str(ir
->init_step
+ir
->nsteps
,sbuf
),
1802 (ir
->init_step
+ir
->nsteps
)*ir
->delta_t
,
1803 gmx_step_str(ir
->init_step
,sbuf2
),
1804 ir
->init_step
*ir
->delta_t
);
1808 fprintf(stderr
,"%s steps, %8.1f ps.\n",
1809 gmx_step_str(ir
->nsteps
,sbuf
),ir
->nsteps
*ir
->delta_t
);
1812 fprintf(fplog
,"\n");
1815 /* Set and write start time */
1816 runtime_start(runtime
);
1817 print_date_and_time(fplog
,cr
->nodeid
,"Started mdrun",runtime
);
1818 wallcycle_start(wcycle
,ewcRUN
);
1820 fprintf(fplog
,"\n");
1822 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
1824 chkpt_ret
=fcCheckPointParallel( cr
->nodeid
,
1826 if ( chkpt_ret
== 0 )
1827 gmx_fatal( 3,__FILE__
,__LINE__
, "Checkpoint error on step %d\n", 0 );
1831 /***********************************************************
1833 * Loop over MD steps
1835 ************************************************************/
1837 /* if rerunMD then read coordinates and velocities from input trajectory */
1840 if (getenv("GMX_FORCE_UPDATE"))
1842 bForceUpdate
= TRUE
;
1845 bNotLastFrame
= read_first_frame(oenv
,&status
,
1846 opt2fn("-rerun",nfile
,fnm
),
1847 &rerun_fr
,TRX_NEED_X
| TRX_READ_V
);
1848 if (rerun_fr
.natoms
!= top_global
->natoms
)
1851 "Number of atoms in trajectory (%d) does not match the "
1852 "run input file (%d)\n",
1853 rerun_fr
.natoms
,top_global
->natoms
);
1855 if (ir
->ePBC
!= epbcNONE
)
1859 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
);
1861 if (max_cutoff2(ir
->ePBC
,rerun_fr
.box
) < sqr(fr
->rlistlong
))
1863 gmx_fatal(FARGS
,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr
.step
,rerun_fr
.time
);
1866 /* Set the shift vectors.
1867 * Necessary here when have a static box different from the tpr box.
1869 calc_shifts(rerun_fr
.box
,fr
->shift_vec
);
1873 /* loop over MD steps or if rerunMD to end of input trajectory */
1875 /* Skip the first Nose-Hoover integration when we get the state from tpx */
1876 bStateFromTPX
= !opt2bSet("-cpi",nfile
,fnm
);
1877 bInitStep
= bFirstStep
&& (bStateFromTPX
|| bVV
);
1878 bStartingFromCpt
= (Flags
& MD_STARTFROMCPT
) && bInitStep
;
1880 bSumEkinhOld
= FALSE
;
1883 step
= ir
->init_step
;
1886 if (ir
->nstlist
== -1)
1888 init_nlistheuristics(&nlh
,bGStatEveryStep
,step
);
1891 bLastStep
= (bRerunMD
|| step_rel
> ir
->nsteps
);
1892 while (!bLastStep
|| (bRerunMD
&& bNotLastFrame
)) {
1894 wallcycle_start(wcycle
,ewcSTEP
);
1896 GMX_MPE_LOG(ev_timestep1
);
1899 if (rerun_fr
.bStep
) {
1900 step
= rerun_fr
.step
;
1901 step_rel
= step
- ir
->init_step
;
1903 if (rerun_fr
.bTime
) {
1913 bLastStep
= (step_rel
== ir
->nsteps
);
1914 t
= t0
+ step
*ir
->delta_t
;
1917 if (ir
->efep
!= efepNO
)
1919 if (bRerunMD
&& rerun_fr
.bLambda
&& (ir
->delta_lambda
!=0))
1921 state_global
->lambda
= rerun_fr
.lambda
;
1925 state_global
->lambda
= lam0
+ step
*ir
->delta_lambda
;
1927 state
->lambda
= state_global
->lambda
;
1928 bDoDHDL
= do_per_step(step
,ir
->nstdhdl
);
1933 update_annealing_target_temp(&(ir
->opts
),t
);
1938 if (!(DOMAINDECOMP(cr
) && !MASTER(cr
)))
1940 for(i
=0; i
<state_global
->natoms
; i
++)
1942 copy_rvec(rerun_fr
.x
[i
],state_global
->x
[i
]);
1946 for(i
=0; i
<state_global
->natoms
; i
++)
1948 copy_rvec(rerun_fr
.v
[i
],state_global
->v
[i
]);
1953 for(i
=0; i
<state_global
->natoms
; i
++)
1955 clear_rvec(state_global
->v
[i
]);
1959 fprintf(stderr
,"\nWARNING: Some frames do not contain velocities.\n"
1960 " Ekin, temperature and pressure are incorrect,\n"
1961 " the virial will be incorrect when constraints are present.\n"
1963 bRerunWarnNoV
= FALSE
;
1967 copy_mat(rerun_fr
.box
,state_global
->box
);
1968 copy_mat(state_global
->box
,state
->box
);
1970 if (vsite
&& (Flags
& MD_RERUN_VSITE
))
1972 if (DOMAINDECOMP(cr
))
1974 gmx_fatal(FARGS
,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
1978 /* Following is necessary because the graph may get out of sync
1979 * with the coordinates if we only have every N'th coordinate set
1981 mk_mshift(fplog
,graph
,fr
->ePBC
,state
->box
,state
->x
);
1982 shift_self(graph
,state
->box
,state
->x
);
1984 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,state
->v
,
1985 top
->idef
.iparams
,top
->idef
.il
,
1986 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
1989 unshift_self(graph
,state
->box
,state
->x
);
1994 /* Stop Center of Mass motion */
1995 bStopCM
= (ir
->comm_mode
!= ecmNO
&& do_per_step(step
,ir
->nstcomm
));
1997 /* Copy back starting coordinates in case we're doing a forcefield scan */
2000 for(ii
=0; (ii
<state
->natoms
); ii
++)
2002 copy_rvec(xcopy
[ii
],state
->x
[ii
]);
2003 copy_rvec(vcopy
[ii
],state
->v
[ii
]);
2005 copy_mat(boxcopy
,state
->box
);
2010 /* for rerun MD always do Neighbour Searching */
2011 bNS
= (bFirstStep
|| ir
->nstlist
!= 0);
2016 /* Determine whether or not to do Neighbour Searching and LR */
2017 bNStList
= (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0);
2019 bNS
= (bFirstStep
|| bExchanged
|| bNStList
||
2020 (ir
->nstlist
== -1 && nlh
.nabnsb
> 0));
2022 if (bNS
&& ir
->nstlist
== -1)
2024 set_nlistheuristics(&nlh
,bFirstStep
|| bExchanged
,step
);
2028 if (terminate_now
> 0 || (terminate_now
< 0 && bNS
))
2033 /* Determine whether or not to update the Born radii if doing GB */
2034 bBornRadii
=bFirstStep
;
2035 if (ir
->implicit_solvent
&& (step
% ir
->nstgbradii
==0))
2040 do_log
= do_per_step(step
,ir
->nstlog
) || bFirstStep
|| bLastStep
;
2041 do_verbose
= bVerbose
&&
2042 (step
% stepout
== 0 || bFirstStep
|| bLastStep
);
2044 if (bNS
&& !(bFirstStep
&& ir
->bContinuation
&& !bRerunMD
))
2048 bMasterState
= TRUE
;
2052 bMasterState
= FALSE
;
2053 /* Correct the new box if it is too skewed */
2054 if (DYNAMIC_BOX(*ir
))
2056 if (correct_box(fplog
,step
,state
->box
,graph
))
2058 bMasterState
= TRUE
;
2061 if (DOMAINDECOMP(cr
) && bMasterState
)
2063 dd_collect_state(cr
->dd
,state
,state_global
);
2067 if (DOMAINDECOMP(cr
))
2069 /* Repartition the domain decomposition */
2070 wallcycle_start(wcycle
,ewcDOMDEC
);
2071 dd_partition_system(fplog
,step
,cr
,
2072 bMasterState
,nstglobalcomm
,
2073 state_global
,top_global
,ir
,
2074 state
,&f
,mdatoms
,top
,fr
,
2075 vsite
,shellfc
,constr
,
2076 nrnb
,wcycle
,do_verbose
);
2077 wallcycle_stop(wcycle
,ewcDOMDEC
);
2081 if (MASTER(cr
) && do_log
&& !bFFscan
)
2083 print_ebin_header(fplog
,step
,t
,state
->lambda
);
2086 if (ir
->efep
!= efepNO
)
2088 update_mdatoms(mdatoms
,state
->lambda
);
2091 if (bRerunMD
&& rerun_fr
.bV
)
2094 /* We need the kinetic energy at minus the half step for determining
2095 * the full step kinetic energy and possibly for T-coupling.*/
2096 /* This may not be quite working correctly yet . . . . */
2097 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2098 wcycle
,enerd
,NULL
,NULL
,NULL
,NULL
,mu_tot
,
2099 constr
,&chkpt
,&terminate
,&terminate_now
,&(nlh
.nabnsb
),state
->box
,
2100 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2101 (CGLO_RERUNMD
| CGLO_GSTAT
| CGLO_TEMPERATURE
) |
2102 (bEkinAveVel
?CGLO_EKINAVEVEL
:0));
2104 clear_mat(force_vir
);
2106 /* Ionize the atoms if necessary */
2109 ionize(fplog
,oenv
,mdatoms
,top_global
,t
,ir
,state
->x
,state
->v
,
2110 mdatoms
->start
,mdatoms
->start
+mdatoms
->homenr
,state
->box
,cr
);
2113 /* Update force field in ffscan program */
2116 if (update_forcefield(fplog
,
2118 mdatoms
->nr
,state
->x
,state
->box
)) {
2119 if (gmx_parallel_env())
2127 GMX_MPE_LOG(ev_timestep2
);
2129 if ((bNS
|| bLastStep
) && (step
> ir
->init_step
) && !bRerunMD
)
2131 bCPT
= (chkpt
> 0 || (bLastStep
&& (Flags
& MD_CONFOUT
)));
2142 /* Determine the pressure:
2143 * always when we want exact averages in the energy file,
2144 * at ns steps when we have pressure coupling,
2145 * otherwise only at energy output steps (set below).
2147 bNstEner
= (bGStatEveryStep
|| do_per_step(step
,ir
->nstcalcenergy
));
2148 bCalcEner
= bNstEner
;
2149 bCalcPres
= (bGStatEveryStep
|| (ir
->epc
!= epcNO
&& bNS
));
2151 /* Do we need global communication ? */
2152 bGStat
= (bCalcEner
|| bStopCM
||
2153 (ir
->nstlist
== -1 && !bRerunMD
&& step
>= nlh
.step_nscheck
));
2155 do_ene
= (do_per_step(step
,ir
->nstenergy
) || bLastStep
);
2157 if (do_ene
|| do_log
)
2164 /* these CGLO_ options remain the same throughout the iteration */
2165 cglo_flags
= ((bRerunMD
? CGLO_RERUNMD
: 0) |
2166 (bEkinAveVel
? CGLO_EKINAVEVEL
: 0) |
2167 (bStopCM
? CGLO_STOPCM
: 0) |
2168 (bGStat
? CGLO_GSTAT
: 0) |
2169 (bNEMD
? CGLO_NEMD
: 0));
2171 force_flags
= (GMX_FORCE_STATECHANGED
|
2172 ((DYNAMIC_BOX(*ir
) || bRerunMD
) ? GMX_FORCE_DYNAMICBOX
: 0) |
2173 GMX_FORCE_ALLFORCES
|
2174 (bNStList
? GMX_FORCE_DOLR
: 0) |
2176 (bCalcEner
? GMX_FORCE_VIRIAL
: 0) |
2177 (bDoDHDL
? GMX_FORCE_DHDL
: 0));
2181 /* Now is the time to relax the shells */
2182 count
=relax_shell_flexcon(fplog
,cr
,bVerbose
,bFFscan
? step
+1 : step
,
2184 bStopCM
,top
,top_global
,
2186 state
,f
,force_vir
,mdatoms
,
2187 nrnb
,wcycle
,graph
,groups
,
2188 shellfc
,fr
,bBornRadii
,t
,mu_tot
,
2189 state
->natoms
,&bConverged
,vsite
,
2200 /* The coordinates (x) are shifted (to get whole molecules)
2202 * This is parallellized as well, and does communication too.
2203 * Check comments in sim_util.c
2206 do_force(fplog
,cr
,ir
,step
,nrnb
,wcycle
,top
,top_global
,groups
,
2207 state
->box
,state
->x
,&state
->hist
,
2208 f
,force_vir
,mdatoms
,enerd
,fcd
,
2209 state
->lambda
,graph
,
2210 fr
,vsite
,mu_tot
,t
,fp_field
,ed
,bBornRadii
,
2211 (bNS
? GMX_FORCE_NS
: 0) | force_flags
);
2214 GMX_BARRIER(cr
->mpi_comm_mygroup
);
2218 mu_aver
= calc_mu_aver(cr
,state
->x
,mdatoms
->chargeA
,
2219 mu_tot
,&top_global
->mols
,mdatoms
,gnx
,grpindex
);
2222 if (bTCR
&& bFirstStep
)
2224 tcr
=init_coupling(fplog
,nfile
,fnm
,cr
,fr
,mdatoms
,&(top
->idef
));
2225 fprintf(fplog
,"Done init_coupling\n");
2229 /* ############### START FIRST UPDATE HALF-STEP ############### */
2231 if (!bStartingFromCpt
&& !bRerunMD
) {
2232 if (bVVAveVel
&& bInitStep
) {
2233 /* if using velocity verlet with full time step Ekin, take the first half step only to compute the
2234 virial for the first step. From there, revert back to the initial coordinates
2235 so that the input is actually the initial step */
2236 copy_rvecn(state
->v
,cbuf
,0,state
->natoms
); /* should make this better for parallelizing? */
2239 /* this is for NHC in the Ekin(t+dt/2) version of vv */
2240 if (!bInitStep
|| !bEkinAveVel
)
2242 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[1],bEkinAveVel
);
2245 update_coords(fplog
,step
,ir
,mdatoms
,state
,
2246 f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2247 ekind
,M
,wcycle
,upd
,bInitStep
,etrtVELOCITY
,cr
,nrnb
,constr
,&top
->idef
);
2249 bIterateFirstHalf
= bIterate
&& (!bInitStep
);
2250 bFirstIterate
= TRUE
;
2251 /* for iterations, we save these vectors, as we will be self-consistently iterating
2253 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
2255 /* save the state */
2256 if (bIterateFirstHalf
) {
2257 copy_coupling_state(state
,bufstate
,ekind
,ekind_save
);
2260 while (bIterateFirstHalf
|| bFirstIterate
)
2262 if (bIterateFirstHalf
)
2264 copy_coupling_state(bufstate
,state
,ekind_save
,ekind
);
2266 if (bIterateFirstHalf
)
2268 if (bFirstIterate
&& bTrotter
)
2270 /* The first time, we need a decent first estimate
2271 of veta(t+dt) to compute the constraints. Do
2272 this by computing the box volume part of the
2273 trotter integration at this time. Nothing else
2274 should be changed by this routine here. If
2275 !(first time), we start with the previous value
2278 veta_save
= state
->veta
;
2279 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[0],FALSE
);
2280 vetanew
= state
->veta
;
2281 state
->veta
= veta_save
;
2286 if ( !bRerunMD
|| rerun_fr
.bV
|| bForceUpdate
) { /* Why is rerun_fr.bV here? Unclear. */
2287 wallcycle_start(wcycle
,ewcUPDATE
);
2290 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
2291 &top
->idef
,shake_vir
,NULL
,
2292 cr
,nrnb
,wcycle
,upd
,constr
,
2293 bInitStep
,TRUE
,bCalcPres
,vetanew
);
2295 if (!bOK
&& !bFFscan
)
2297 gmx_fatal(FARGS
,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2302 { /* Need to unshift here if a do_force has been
2303 called in the previous step */
2304 unshift_self(graph
,state
->box
,state
->x
);
2307 /* if bEkinAveVel, then compute ekin and shake_vir, otherwise, just compute shake_vir */
2308 /* first step, we need to compute the virial */
2310 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2311 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
2312 constr
,&chkpt
,&terminate
,&terminate_now
,&(nlh
.nabnsb
),state
->box
,
2313 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2314 (cglo_flags
| CGLO_ENERGY
) |
2315 ((bIterateFirstHalf
| bInitStep
| IR_NPT_TROTTER(ir
)) ? (cglo_flags
| CGLO_PRESSURE
):0) |
2316 (((bEkinAveVel
&&(!bInitStep
)) || (!bEkinAveVel
&& IR_NPT_TROTTER(ir
)))? (cglo_flags
| CGLO_TEMPERATURE
):0) |
2317 ((bEkinAveVel
|| (!bEkinAveVel
&& IR_NPT_TROTTER(ir
))) ? CGLO_EKINAVEVEL
: 0) |
2318 (bIterate
? CGLO_ITERATE
: 0) |
2319 (bFirstIterate
? CGLO_FIRSTITERATE
: 0) |
2320 (cglo_flags
| CGLO_SCALEEKIN
)
2323 /* explanation of above:
2324 a) We compute Ekin at the full time step
2325 if 1) we are using the AveVel Ekin, and it's not the
2326 initial step, or 2) if we are using AveEkin, but need the full
2327 time step kinetic energy for the pressure.
2328 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
2329 EkinAveVel because it's needed for the pressure */
2331 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
2334 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[2],bEkinAveVel
);
2337 if (done_iterating(cr
,fplog
,&bFirstIterate
,&bIterateFirstHalf
,state
->veta
,&vetanew
,0)) break;
2340 if (bTrotter
&& !bInitStep
) {
2341 copy_mat(shake_vir
,state
->vir_prev
);
2342 if (IR_NVT_TROTTER(ir
) && bEkinAveVel
) {
2343 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
2344 enerd
->term
[F_TEMP
] = sum_ekin(&(ir
->opts
),ekind
,NULL
,bEkinAveVel
,FALSE
,FALSE
);
2345 enerd
->term
[F_EKIN
] = trace(ekind
->ekin
);
2348 /* if it's the initial step, we performed this first step just to get the constraint virial */
2349 if (bInitStep
&& bVVAveVel
) {
2350 copy_rvecn(cbuf
,state
->v
,0,state
->natoms
);
2353 if (fr
->bSepDVDL
&& fplog
&& do_log
)
2355 fprintf(fplog
,sepdvdlformat
,"Constraint",0.0,dvdl
);
2357 enerd
->term
[F_DHDL_CON
] += dvdl
;
2359 wallcycle_stop(wcycle
,ewcUPDATE
);
2360 GMX_MPE_LOG(ev_timestep1
);
2362 /* MRS -- done iterating -- compute the conserved quantity */
2366 if (IR_NVT_TROTTER(ir
) || IR_NPT_TROTTER(ir
))
2369 NPT_energy(ir
,state
->nosehoover_xi
,state
->nosehoover_vxi
,state
->veta
, state
->box
, &MassQ
);
2370 if ((ir
->eDispCorr
!= edispcEnerPres
) && (ir
->eDispCorr
!= edispcAllEnerPres
))
2372 last_conserved
-= enerd
->term
[F_DISPCORR
];
2376 last_ekin
= enerd
->term
[F_EKIN
]; /* does this get preserved through checkpointing? */
2380 /* ######## END FIRST UPDATE STEP ############## */
2381 /* ######## If doing VV, we now have v(dt) ###### */
2383 /* ################## START TRAJECTORY OUTPUT ################# */
2385 /* Now we have the energies and forces corresponding to the
2386 * coordinates at time t. We must output all of this before
2388 * for RerunMD t is read from input trajectory
2390 GMX_MPE_LOG(ev_output_start
);
2392 bX
= do_per_step(step
,ir
->nstxout
);
2393 bV
= do_per_step(step
,ir
->nstvout
);
2394 bF
= do_per_step(step
,ir
->nstfout
);
2395 bXTC
= do_per_step(step
,ir
->nstxtcout
);
2399 fcReportProgress( ir
->nsteps
, step
);
2401 bX
= bX
|| bLastStep
; /*enforce writing positions and velocities
2403 bV
= bV
|| bLastStep
;
2405 int nthreads
=(cr
->nthreads
==0 ? 1 : cr
->nthreads
);
2406 int nnodes
=(cr
->nnodes
==0 ? 1 : cr
->nnodes
);
2409 /*Gromacs drives checkpointing; no ||
2410 fcCheckPointPendingThreads(cr->nodeid,
2412 /* sync bCPT and fc record-keeping */
2413 if (bCPT
&& MASTER(cr
))
2414 fcRequestCheckPoint();
2419 if (bX
|| bV
|| bF
|| bXTC
|| bCPT
)
2421 wallcycle_start(wcycle
,ewcTRAJ
);
2424 if (state
->flags
& (1<<estLD_RNG
))
2426 get_stochd_state(upd
,state
);
2432 state_global
->ekinstate
.bUpToDate
= FALSE
;
2436 update_ekinstate(&state_global
->ekinstate
,ekind
);
2437 state_global
->ekinstate
.bUpToDate
= TRUE
;
2439 update_energyhistory(&state_global
->enerhist
,mdebin
);
2442 write_traj(fplog
,cr
,fp_trn
,bX
,bV
,bF
,fp_xtc
,bXTC
,ir
->xtcprec
,
2443 fn_cpt
,bCPT
,top_global
,ir
->eI
,ir
->simulation_part
,
2444 step
,t
,state
,state_global
,f
,f_global
,&n_xtc
,&x_xtc
);
2446 if (bLastStep
&& step_rel
== ir
->nsteps
&&
2447 (Flags
& MD_CONFOUT
) && MASTER(cr
) &&
2448 !bRerunMD
&& !bFFscan
)
2450 /* x and v have been collected in write_traj */
2451 fprintf(stderr
,"\nWriting final coordinates.\n");
2452 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
&&
2455 /* Make molecules whole only for confout writing */
2456 do_pbc_mtop(fplog
,ir
->ePBC
,state
->box
,top_global
,state_global
->x
);
2458 write_sto_conf_mtop(ftp2fn(efSTO
,nfile
,fnm
),
2459 *top_global
->name
,top_global
,
2460 state_global
->x
,state_global
->v
,
2461 ir
->ePBC
,state
->box
);
2464 wallcycle_stop(wcycle
,ewcTRAJ
);
2466 GMX_MPE_LOG(ev_output_finish
);
2468 /* kludge -- virial is lost with restart for NPT control. Must restart */
2469 if (bStartingFromCpt
&& bVV
)
2471 copy_mat(state
->vir_prev
,shake_vir
);
2473 /* ################## END TRAJECTORY OUTPUT ################ */
2475 /* Determine the pressure:
2476 * always when we want exact averages in the energy file,
2477 * at ns steps when we have pressure coupling,
2478 * otherwise only at energy output steps (set below).
2481 bNstEner
= (bGStatEveryStep
|| do_per_step(step
,ir
->nstcalcenergy
));
2482 bCalcEner
= bNstEner
;
2483 bCalcPres
= (bGStatEveryStep
|| (ir
->epc
!= epcNO
&& bNS
));
2485 /* Do we need global communication ? */
2486 bGStat
= (bGStatEveryStep
|| bStopCM
|| bNS
||
2487 (ir
->nstlist
== -1 && !bRerunMD
&& step
>= nlh
.step_nscheck
));
2489 do_ene
= (do_per_step(step
,ir
->nstenergy
) || bLastStep
);
2491 if (do_ene
|| do_log
)
2497 bFirstIterate
= TRUE
;
2499 /* for iterations, we save these vectors, as we will be redoing the calculations */
2502 copy_coupling_state(state
,bufstate
,ekind
,ekind_save
);
2505 while (bIterate
|| bFirstIterate
)
2507 /* We now restore these vectors to redo the calculation with improved extended variables */
2510 copy_coupling_state(bufstate
,state
,ekind_save
,ekind
);
2513 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
2514 so scroll down for that logic */
2516 /* ######### START SECOND UPDATE STEP ################# */
2518 GMX_MPE_LOG(ev_update_start
);
2520 if (!bRerunMD
|| rerun_fr
.bV
|| bForceUpdate
)
2522 wallcycle_start(wcycle
,ewcUPDATE
);
2526 /* Box is changed in update() when we do pressure coupling,
2527 * but we should still use the old box for energy corrections and when
2528 * writing it to the energy file, so it matches the trajectory files for
2529 * the same timestep above. Make a copy in a separate array.
2532 copy_mat(state
->box
,lastbox
);
2534 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS*/
2539 tracevir
= trace(shake_vir
);
2543 /* we use a new value of scalevir to converge the iterations faster */
2544 scalevir
= tracevir
/trace(shake_vir
);
2545 msmul(shake_vir
,scalevir
,shake_vir
);
2546 m_add(force_vir
,shake_vir
,total_vir
);
2547 clear_mat(shake_vir
);
2549 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[3],bEkinAveVel
);
2551 /* We can only do Berendsen coupling after we have summed the kinetic
2552 * energy or virial. Since the happens in global_state after update,
2553 * we should only do it at step % nstlist = 1 with bGStatEveryStep=FALSE.
2556 update_extended(fplog
,step
,ir
,mdatoms
,state
,ekind
,pcoupl_mu
,M
,wcycle
,
2557 upd
,bInitStep
,FALSE
,&MassQ
);
2559 /* velocity (for VV) */
2560 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2561 ekind
,M
,wcycle
,upd
,FALSE
,etrtVELOCITY
,cr
,nrnb
,constr
,&top
->idef
);
2563 /* above, initialize just copies ekinh into ekin, it doesn't copy
2564 /* position (for VV), and entire integrator for MD */
2568 copy_rvecn(state
->x
,cbuf
,0,state
->natoms
);
2571 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2572 ekind
,M
,wcycle
,upd
,bInitStep
,etrtPOSITION
,cr
,nrnb
,constr
,&top
->idef
);
2574 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
2575 &top
->idef
,shake_vir
,force_vir
,
2576 cr
,nrnb
,wcycle
,upd
,constr
,
2577 bInitStep
,FALSE
,bCalcPres
,state
->veta
);
2581 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2582 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
2583 constr
,&chkpt
,&terminate
,&terminate_now
,&(nlh
.nabnsb
),lastbox
,
2584 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2585 (cglo_flags
| CGLO_TEMPERATURE
) |
2586 (cglo_flags
& ~CGLO_PRESSURE
) |
2587 (cglo_flags
& ~CGLO_ENERGY
)
2589 trotter_update(ir
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
[4],bEkinAveVel
);
2590 copy_rvecn(cbuf
,state
->x
,0,state
->natoms
);
2591 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2592 ekind
,M
,wcycle
,upd
,bInitStep
,etrtPOSITION
,cr
,nrnb
,constr
,&top
->idef
);
2593 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
2594 /* are the small terms in the shake_vir here due
2595 * to numerical errors, or are they important
2596 * physically? I'm thinking they are just errors, but not completely sure. */
2597 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
2598 &top
->idef
,tmp_vir
,force_vir
,
2599 cr
,nrnb
,wcycle
,upd
,constr
,
2600 bInitStep
,FALSE
,bCalcPres
,state
->veta
);
2601 m_add(shake_vir
,tmp_vir
,shake_vir
);
2604 if (!bOK
&& !bFFscan
)
2606 gmx_fatal(FARGS
,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2609 if (fr
->bSepDVDL
&& fplog
&& do_log
)
2611 fprintf(fplog
,sepdvdlformat
,"Constraint",0.0,dvdl
);
2613 enerd
->term
[F_DHDL_CON
] += dvdl
;
2614 wallcycle_stop(wcycle
,ewcUPDATE
);
2618 /* Need to unshift here */
2619 unshift_self(graph
,state
->box
,state
->x
);
2622 GMX_BARRIER(cr
->mpi_comm_mygroup
);
2623 GMX_MPE_LOG(ev_update_finish
);
2627 wallcycle_start(wcycle
,ewcVSITECONSTR
);
2630 shift_self(graph
,state
->box
,state
->x
);
2632 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,state
->v
,
2633 top
->idef
.iparams
,top
->idef
.il
,
2634 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
2638 unshift_self(graph
,state
->box
,state
->x
);
2640 wallcycle_stop(wcycle
,ewcVSITECONSTR
);
2643 /* ############## IF NOT VV, CALCULATE EKIN HERE - ALWAYS CALCULATE PRESSURE ############ */
2644 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2645 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
2646 constr
,&chkpt
,&terminate
,&terminate_now
,&(nlh
.nabnsb
),lastbox
,
2647 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2648 (!bVV
? (cglo_flags
| CGLO_ENERGY
):0) |
2649 (cglo_flags
| CGLO_PRESSURE
) |
2650 (!bVV
? (cglo_flags
| CGLO_TEMPERATURE
):0) |
2651 (bIterate
? CGLO_ITERATE
: 0) |
2652 (bFirstIterate
? CGLO_FIRSTITERATE
: 0)
2654 /* bIterate is set to keep it from eliminating the old ekinh */
2655 /* ############# END CALC EKIN AND PRESSURE ################# */
2657 if (done_iterating(cr
,fplog
,&bFirstIterate
,&bIterate
,trace(shake_vir
),&tracevir
,1)) break;
2660 update_box(fplog
,step
,ir
,mdatoms
,state
,graph
,f
,
2661 ir
->nstlist
==-1 ? &nlh
.scale_tot
: NULL
,pcoupl_mu
,nrnb
,wcycle
,upd
,bInitStep
,FALSE
);
2663 /* ################# END UPDATE STEP 2 ################# */
2664 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
2666 /* Determine the wallclock run time up till now */
2667 run_time
= (double)time(NULL
) - (double)runtime
->real
;
2669 /* Check whether everything is still allright */
2670 if ((bGotTermSignal
|| bGotUsr1Signal
) && !bHandledSignal
)
2672 if (bGotTermSignal
|| ir
->nstlist
== 0)
2682 terminate_now
= terminate
;
2687 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2688 bGotTermSignal
? "TERM" : "USR1",
2689 terminate
==-1 ? "NS " : "");
2693 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2694 bGotTermSignal
? "TERM" : "USR1",
2695 terminate
==-1 ? "NS " : "");
2697 bHandledSignal
=TRUE
;
2699 else if (MASTER(cr
) && (bNS
|| ir
->nstlist
<= 0) &&
2700 (max_hours
> 0 && run_time
> max_hours
*60.0*60.0*0.99) &&
2703 /* Signal to terminate the run */
2704 terminate
= (ir
->nstlist
== 0 ? 1 : -1);
2707 terminate_now
= terminate
;
2711 fprintf(fplog
,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step
,sbuf
),max_hours
*0.99);
2713 fprintf(stderr
, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step
,sbuf
),max_hours
*0.99);
2716 if (ir
->nstlist
== -1 && !bRerunMD
)
2718 /* When bGStatEveryStep=FALSE, global_stat is only called
2719 * when we check the atom displacements, not at NS steps.
2720 * This means that also the bonded interaction count check is not
2721 * performed immediately after NS. Therefore a few MD steps could
2722 * be performed with missing interactions.
2723 * But wrong energies are never written to file,
2724 * since energies are only written after global_stat
2727 if (step
>= nlh
.step_nscheck
)
2729 nlh
.nabnsb
= natoms_beyond_ns_buffer(ir
,fr
,&top
->cgs
,
2730 nlh
.scale_tot
,state
->x
);
2734 /* This is not necessarily true,
2735 * but step_nscheck is determined quite conservatively.
2741 /* In parallel we only have to check for checkpointing in steps
2742 * where we do global communication,
2743 * otherwise the other nodes don't know.
2745 if (MASTER(cr
) && ((bGStat
|| !PAR(cr
)) &&
2748 run_time
>= nchkpt
*cpt_period
*60.0)))
2757 /* The coordinates (x) were unshifted in update */
2758 if (bFFscan
&& (shellfc
==NULL
|| bConverged
))
2760 if (print_forcefield(fplog
,enerd
->term
,mdatoms
->homenr
,
2762 &(top_global
->mols
),mdatoms
->massT
,pres
))
2764 if (gmx_parallel_env()) {
2767 fprintf(stderr
,"\n");
2773 /* We will not sum ekinh_old,
2774 * so signal that we still have to do it.
2776 bSumEkinhOld
= TRUE
;
2781 /* Only do GCT when the relaxation of shells (minimization) has converged,
2782 * otherwise we might be coupling to bogus energies.
2783 * In parallel we must always do this, because the other sims might
2787 /* Since this is called with the new coordinates state->x, I assume
2788 * we want the new box state->box too. / EL 20040121
2790 do_coupling(fplog
,oenv
,nfile
,fnm
,tcr
,t
,step
,enerd
->term
,fr
,
2792 mdatoms
,&(top
->idef
),mu_aver
,
2793 top_global
->mols
.nr
,cr
,
2794 state
->box
,total_vir
,pres
,
2795 mu_tot
,state
->x
,f
,bConverged
);
2799 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
2801 sum_dhdl(enerd
,state
->lambda
,ir
);
2802 /* use the directly determined last velocity, not actually the averaged half steps */
2803 if (bTrotter
&& bEkinAveVel
) {
2804 enerd
->term
[F_EKIN
] = last_ekin
;
2806 enerd
->term
[F_ETOT
] = enerd
->term
[F_EPOT
] + enerd
->term
[F_EKIN
];
2815 if (IR_NVT_TROTTER(ir
)) {
2816 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] + last_conserved
;
2818 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] +
2819 NPT_energy(ir
,state
->nosehoover_xi
,
2820 state
->nosehoover_vxi
,state
->veta
,state
->box
,&MassQ
);
2824 enerd
->term
[F_ECONSERVED
] =
2825 enerd
->term
[F_ETOT
] + vrescale_energy(&(ir
->opts
),
2826 state
->therm_integral
);
2832 /* Check for excessively large energies */
2836 real etot_max
= 1e200
;
2838 real etot_max
= 1e30
;
2840 if (fabs(enerd
->term
[F_ETOT
]) > etot_max
)
2842 fprintf(stderr
,"Energy too large (%g), giving up\n",
2843 enerd
->term
[F_ETOT
]);
2846 /* ######### END PREPARING EDR OUTPUT ########### */
2848 /* Time for performance */
2849 if (((step
% stepout
) == 0) || bLastStep
)
2851 runtime_upd_proc(runtime
);
2861 upd_mdebin(mdebin
,bDoDHDL
? fp_dhdl
: NULL
,TRUE
,
2862 t
,mdatoms
->tmass
,enerd
,state
,lastbox
,
2863 shake_vir
,force_vir
,total_vir
,pres
,
2864 ekind
,mu_tot
,constr
);
2868 upd_mdebin_step(mdebin
);
2871 do_dr
= do_per_step(step
,ir
->nstdisreout
);
2872 do_or
= do_per_step(step
,ir
->nstorireout
);
2874 print_ebin(fp_ene
,do_ene
,do_dr
,do_or
,do_log
?fplog
:NULL
,step
,t
,
2875 eprNORMAL
,bCompact
,mdebin
,fcd
,groups
,&(ir
->opts
));
2877 if (ir
->ePull
!= epullNO
)
2879 pull_print_output(ir
->pull
,step
,t
);
2882 if (do_per_step(step
,ir
->nstlog
))
2884 if(fflush(fplog
) != 0)
2886 gmx_fatal(FARGS
,"Cannot flush logfile - maybe you are out of quota?");
2892 /* Remaining runtime */
2893 if (MULTIMASTER(cr
) && do_verbose
)
2897 fprintf(stderr
,"\n");
2899 print_time(stderr
,runtime
,step
,ir
);
2902 /* Replica exchange */
2904 if ((repl_ex_nst
> 0) && (step
> 0) && !bLastStep
&&
2905 do_per_step(step
,repl_ex_nst
))
2907 bExchanged
= replica_exchange(fplog
,cr
,repl_ex
,
2908 state_global
,enerd
->term
,
2911 if (bExchanged
&& PAR(cr
))
2913 if (DOMAINDECOMP(cr
))
2915 dd_partition_system(fplog
,step
,cr
,TRUE
,1,
2916 state_global
,top_global
,ir
,
2917 state
,&f
,mdatoms
,top
,fr
,
2918 vsite
,shellfc
,constr
,
2923 bcast_state(cr
,state
,FALSE
);
2929 bStartingFromCpt
= FALSE
;
2931 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
2932 /* Complicated conditional when bGStatEveryStep=FALSE.
2933 * We can not just use bGStat, since then the simulation results
2934 * would depend on nstenergy and nstlog or step_nscheck.
2936 if (((state
->flags
& (1<<estPRES_PREV
)) || (state
->flags
& (1<<estVIR_PREV
))) &&
2938 (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0) ||
2939 (ir
->nstlist
< 0 && nlh
.nabnsb
> 0) ||
2940 (ir
->nstlist
== 0 && bGStat
)))
2942 /* Store the pressure in t_state for pressure coupling
2943 * at the next MD step.
2945 if (state
->flags
& (1<<estPRES_PREV
))
2947 copy_mat(pres
,state
->pres_prev
);
2951 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
2955 /* read next frame from input trajectory */
2956 bNotLastFrame
= read_next_frame(oenv
,status
,&rerun_fr
);
2959 if (!bRerunMD
|| !rerun_fr
.bStep
)
2961 /* increase the MD step number */
2966 cycles
= wallcycle_stop(wcycle
,ewcSTEP
);
2967 if (DOMAINDECOMP(cr
) && wcycle
)
2969 dd_cycles_add(cr
->dd
,cycles
,ddCyclStep
);
2972 if (step_rel
== wcycle_get_reset_counters(wcycle
))
2974 /* Reset all the counters related to performance over the run */
2975 reset_all_counters(fplog
,cr
,step
,&step_rel
,ir
,wcycle
,nrnb
,runtime
);
2976 wcycle_set_reset_counters(wcycle
, 0);
2979 /* End of main MD loop */
2983 runtime_end(runtime
);
2990 if (!(cr
->duty
& DUTY_PME
))
2992 /* Tell the PME only node to finish */
2998 if (bGStatEveryStep
&& !bRerunMD
)
3000 print_ebin(fp_ene
,FALSE
,FALSE
,FALSE
,fplog
,step
,t
,
3001 eprAVER
,FALSE
,mdebin
,fcd
,groups
,&(ir
->opts
));
3011 gmx_fio_fclose(fp_dhdl
);
3015 gmx_fio_fclose(fp_field
);
3020 if (ir
->nstlist
== -1 && nlh
.nns
> 0 && fplog
)
3022 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
)));
3023 fprintf(fplog
,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh
.ab
/nlh
.nns
);
3026 if (shellfc
&& fplog
)
3028 fprintf(fplog
,"Fraction of iterations that converged: %.2f %%\n",
3029 (nconverged
*100.0)/step_rel
);
3030 fprintf(fplog
,"Average number of force evaluations per MD step: %.2f\n\n",
3034 if (repl_ex_nst
> 0 && MASTER(cr
))
3036 print_replica_exchange_statistics(fplog
,repl_ex
);
3039 runtime
->nsteps_done
= step_rel
;