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
61 #include "pull_rotation.h"
75 #include "mpelogging.h"
82 #include "compute_io.h"
84 #include "checkpoint.h"
85 #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
} };
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
;
133 struct mdrunner_arglist
147 const char *dddlb_opt
;
160 int ret
; /* return value */
164 static void mdrunner_start_fn(void *arg
)
166 struct mdrunner_arglist
*mda
=(struct mdrunner_arglist
*)arg
;
167 struct mdrunner_arglist mc
=*mda
; /* copy the arg list to make sure
168 that it's thread-local. This doesn't
169 copy pointed-to items, of course,
170 but those are all const. */
171 t_commrec
*cr
; /* we need a local version of this */
173 t_filenm
*fnm
=dup_tfn(mc
.nfile
, mc
.fnm
);
175 cr
=init_par_threads(mc
.cr
);
182 mda
->ret
=mdrunner(fplog
, cr
, mc
.nfile
, mc
.fnm
, mc
.oenv
, mc
.bVerbose
,
183 mc
.bCompact
, mc
.nstglobalcomm
,
184 mc
.ddxyz
, mc
.dd_node_order
, mc
.rdd
,
185 mc
.rconstr
, mc
.dddlb_opt
, mc
.dlb_scale
,
186 mc
.ddcsx
, mc
.ddcsy
, mc
.ddcsz
, mc
.nstepout
, mc
.nmultisim
,
187 mc
.repl_ex_nst
, mc
.repl_ex_seed
, mc
.pforce
,
188 mc
.cpt_period
, mc
.max_hours
, mc
.Flags
);
193 int mdrunner_threads(int nthreads
,
194 FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
195 const output_env_t oenv
, bool bVerbose
,bool bCompact
,
197 ivec ddxyz
,int dd_node_order
,real rdd
,real rconstr
,
198 const char *dddlb_opt
,real dlb_scale
,
199 const char *ddcsx
,const char *ddcsy
,const char *ddcsz
,
200 int nstepout
,int nmultisim
,int repl_ex_nst
,
201 int repl_ex_seed
, real pforce
,real cpt_period
,
202 real max_hours
, unsigned long Flags
)
205 /* first check whether we even need to start tMPI */
208 ret
=mdrunner(fplog
, cr
, nfile
, fnm
, oenv
, bVerbose
, bCompact
,
210 ddxyz
, dd_node_order
, rdd
, rconstr
, dddlb_opt
, dlb_scale
,
211 ddcsx
, ddcsy
, ddcsz
, nstepout
, nmultisim
, repl_ex_nst
,
212 repl_ex_seed
, pforce
, cpt_period
, max_hours
, Flags
);
217 struct mdrunner_arglist mda
;
218 /* fill the data structure to pass as void pointer to thread start fn */
224 mda
.bVerbose
=bVerbose
;
225 mda
.bCompact
=bCompact
;
226 mda
.nstglobalcomm
=nstglobalcomm
;
227 mda
.ddxyz
[XX
]=ddxyz
[XX
];
228 mda
.ddxyz
[YY
]=ddxyz
[YY
];
229 mda
.ddxyz
[ZZ
]=ddxyz
[ZZ
];
230 mda
.dd_node_order
=dd_node_order
;
233 mda
.dddlb_opt
=dddlb_opt
;
234 mda
.dlb_scale
=dlb_scale
;
238 mda
.nstepout
=nstepout
;
239 mda
.nmultisim
=nmultisim
;
240 mda
.repl_ex_nst
=repl_ex_nst
;
241 mda
.repl_ex_seed
=repl_ex_seed
;
243 mda
.cpt_period
=cpt_period
;
244 mda
.max_hours
=max_hours
;
247 fprintf(stderr
, "Starting %d threads\n",nthreads
);
249 tMPI_Init_fn(nthreads
, mdrunner_start_fn
, (void*)(&mda
) );
253 gmx_comm("Multiple threads requested but not compiled with threads");
260 int mdrunner(FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
261 const output_env_t oenv
, bool bVerbose
,bool bCompact
,
263 ivec ddxyz
,int dd_node_order
,real rdd
,real rconstr
,
264 const char *dddlb_opt
,real dlb_scale
,
265 const char *ddcsx
,const char *ddcsy
,const char *ddcsz
,
266 int nstepout
,int nmultisim
,int repl_ex_nst
,int repl_ex_seed
,
267 real pforce
,real cpt_period
,real max_hours
,
270 double nodetime
=0,realtime
;
271 t_inputrec
*inputrec
;
278 gmx_mtop_t
*mtop
=NULL
;
279 t_mdatoms
*mdatoms
=NULL
;
283 gmx_pme_t
*pmedata
=NULL
;
284 gmx_vsite_t
*vsite
=NULL
;
286 int i
,m
,nChargePerturbed
=-1,status
,nalloc
;
288 gmx_wallcycle_t wcycle
;
289 bool bReadRNG
,bReadEkin
;
291 gmx_runtime_t runtime
;
293 gmx_step_t reset_counters
;
296 /* Essential dynamics */
297 if (opt2bSet("-ei",nfile
,fnm
))
299 /* Open input and output files, allocate space for ED data structure */
300 ed
= ed_open(nfile
,fnm
,cr
);
308 if (bVerbose
&& SIMMASTER(cr
))
310 fprintf(stderr
,"Getting Loaded...\n");
313 if (Flags
& MD_APPENDFILES
)
320 /* The master thread on the master node reads from disk,
321 * then distributes everything to the other processors.
324 list
= (SIMMASTER(cr
) && !(Flags
& MD_APPENDFILES
)) ? (LIST_SCALARS
| LIST_INPUTREC
) : 0;
327 init_parallel(fplog
, opt2fn_master("-s",nfile
,fnm
,cr
),cr
,
328 inputrec
,mtop
,state
,list
);
333 /* Read a file for a single processor */
335 init_single(fplog
,inputrec
,ftp2fn(efTPX
,nfile
,fnm
),mtop
,state
);
337 if (!EEL_PME(inputrec
->coulombtype
) || (Flags
& MD_PARTDEC
))
343 fcRegisterSteps(inputrec
->nsteps
,inputrec
->init_step
);
346 /* NMR restraints must be initialized before load_checkpoint,
347 * since with time averaging the history is added to t_state.
348 * For proper consistency check we therefore need to extend
350 * So the PME-only nodes (if present) will also initialize
351 * the distance restraints.
355 /* This needs to be called before read_checkpoint to extend the state */
356 init_disres(fplog
,mtop
,inputrec
,cr
,Flags
& MD_PARTDEC
,fcd
,state
);
358 if (gmx_mtop_ftype_count(mtop
,F_ORIRES
) > 0)
360 if (PAR(cr
) && !(Flags
& MD_PARTDEC
))
362 gmx_fatal(FARGS
,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
364 /* Orientation restraints */
367 init_orires(fplog
,mtop
,state
->x
,inputrec
,cr
->ms
,&(fcd
->orires
),
372 if (DEFORM(*inputrec
))
374 /* Store the deform reference box before reading the checkpoint */
377 copy_mat(state
->box
,box
);
381 gmx_bcast(sizeof(box
),box
,cr
);
383 /* Because we do not have the update struct available yet
384 * in which the reference values should be stored,
385 * we store them temporarily in static variables.
386 * This should be thread safe, since they are only written once
387 * and with identical values.
390 tMPI_Thread_mutex_lock(&box_mutex
);
392 init_step_tpx
= inputrec
->init_step
;
393 copy_mat(box
,box_tpx
);
395 tMPI_Thread_mutex_unlock(&box_mutex
);
399 if (opt2bSet("-cpi",nfile
,fnm
))
401 /* Check if checkpoint file exists before doing continuation.
402 * This way we can use identical input options for the first and subsequent runs...
404 if( gmx_fexist_master(opt2fn_master("-cpi",nfile
,fnm
,cr
),cr
) )
406 load_checkpoint(opt2fn_master("-cpi",nfile
,fnm
,cr
),fplog
,
407 cr
,Flags
& MD_PARTDEC
,ddxyz
,
408 inputrec
,state
,&bReadRNG
,&bReadEkin
,
409 (Flags
& MD_APPENDFILES
));
413 Flags
|= MD_READ_RNG
;
417 Flags
|= MD_READ_EKIN
;
422 if (MASTER(cr
) && (Flags
& MD_APPENDFILES
))
424 fplog
= gmx_log_open(ftp2fn(efLOG
,nfile
,fnm
),cr
,!(Flags
& MD_SEPPOT
),
430 copy_mat(state
->box
,box
);
435 gmx_bcast(sizeof(box
),box
,cr
);
438 if (bVerbose
&& SIMMASTER(cr
))
440 fprintf(stderr
,"Loaded with Money\n\n");
443 if (PAR(cr
) && !((Flags
& MD_PARTDEC
) || EI_TPI(inputrec
->eI
)))
445 cr
->dd
= init_domain_decomposition(fplog
,cr
,Flags
,ddxyz
,rdd
,rconstr
,
452 make_dd_communicators(fplog
,cr
,dd_node_order
);
454 /* Set overallocation to avoid frequent reallocation of arrays */
455 set_over_alloc_dd(TRUE
);
459 cr
->duty
= (DUTY_PP
| DUTY_PME
);
460 npme_major
= cr
->nnodes
;
462 if (inputrec
->ePBC
== epbcSCREW
)
465 "pbc=%s is only implemented with domain decomposition",
466 epbc_names
[inputrec
->ePBC
]);
472 /* After possible communicator splitting in make_dd_communicators.
473 * we can set up the intra/inter node communication.
475 gmx_setup_nodecomm(fplog
,cr
);
478 wcycle
= wallcycle_init(fplog
,cr
);
481 /* Master synchronizes its value of reset_counters with all nodes
482 * including PME only nodes */
483 reset_counters
= wcycle_get_reset_counters(wcycle
);
484 gmx_bcast_sim(sizeof(reset_counters
),&reset_counters
,cr
);
485 wcycle_set_reset_counters(wcycle
, reset_counters
);
490 if (cr
->duty
& DUTY_PP
)
492 /* For domain decomposition we allocate dynamically
493 * in dd_partition_system.
495 if (DOMAINDECOMP(cr
))
497 bcast_state_setup(cr
,state
);
507 bcast_state(cr
,state
,TRUE
);
511 /* Dihedral Restraints */
512 if (gmx_mtop_ftype_count(mtop
,F_DIHRES
) > 0)
514 init_dihres(fplog
,mtop
,inputrec
,fcd
);
517 /* Initiate forcerecord */
519 init_forcerec(fplog
,oenv
,fr
,fcd
,inputrec
,mtop
,cr
,box
,FALSE
,
520 opt2fn("-table",nfile
,fnm
),
521 opt2fn("-tablep",nfile
,fnm
),
522 opt2fn("-tableb",nfile
,fnm
),FALSE
,pforce
);
524 /* version for PCA_NOT_READ_NODE (see md.c) */
525 /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
526 "nofile","nofile","nofile",FALSE,pforce);
528 fr
->bSepDVDL
= ((Flags
& MD_SEPPOT
) == MD_SEPPOT
);
530 /* Initialize QM-MM */
533 init_QMMMrec(cr
,box
,mtop
,inputrec
,fr
);
536 /* Initialize the mdatoms structure.
537 * mdatoms is not filled with atom data,
538 * as this can not be done now with domain decomposition.
540 mdatoms
= init_mdatoms(fplog
,mtop
,inputrec
->efep
!=efepNO
);
542 /* Initialize the virtual site communication */
543 vsite
= init_vsite(mtop
,cr
);
545 calc_shifts(box
,fr
->shift_vec
);
547 /* With periodic molecules the charge groups should be whole at start up
548 * and the virtual sites should not be far from their proper positions.
550 if (!inputrec
->bContinuation
&& MASTER(cr
) &&
551 !(inputrec
->ePBC
!= epbcNONE
&& inputrec
->bPeriodicMols
))
553 /* Make molecules whole at start of run */
554 if (fr
->ePBC
!= epbcNONE
)
556 do_pbc_first_mtop(fplog
,inputrec
->ePBC
,box
,mtop
,state
->x
);
560 /* Correct initial vsite positions are required
561 * for the initial distribution in the domain decomposition
562 * and for the initial shell prediction.
564 construct_vsites_mtop(fplog
,vsite
,mtop
,state
->x
);
568 /* Initiate PPPM if necessary */
569 if (fr
->eeltype
== eelPPPM
)
571 if (mdatoms
->nChargePerturbed
)
573 gmx_fatal(FARGS
,"Free energy with %s is not implemented",
574 eel_names
[fr
->eeltype
]);
576 status
= gmx_pppm_init(fplog
,cr
,oenv
,FALSE
,TRUE
,box
,
577 getenv("GMXGHAT"),inputrec
, (Flags
& MD_REPRODUCIBLE
));
580 gmx_fatal(FARGS
,"Error %d initializing PPPM",status
);
584 if (EEL_PME(fr
->eeltype
))
586 ewaldcoeff
= fr
->ewaldcoeff
;
587 pmedata
= &fr
->pmedata
;
596 /* This is a PME only node */
598 /* We don't need the state */
601 ewaldcoeff
= calc_ewaldcoeff(inputrec
->rcoulomb
, inputrec
->ewald_rtol
);
605 /* Initiate PME if necessary,
606 * either on all nodes or on dedicated PME nodes only. */
607 if (EEL_PME(inputrec
->coulombtype
))
611 nChargePerturbed
= mdatoms
->nChargePerturbed
;
613 if (cr
->npmenodes
> 0)
615 /* The PME only nodes need to know nChargePerturbed */
616 gmx_bcast_sim(sizeof(nChargePerturbed
),&nChargePerturbed
,cr
);
618 if (cr
->duty
& DUTY_PME
)
620 status
= gmx_pme_init(pmedata
,cr
,npme_major
,inputrec
,
621 mtop
? mtop
->natoms
: 0,nChargePerturbed
,
622 (Flags
& MD_REPRODUCIBLE
));
624 gmx_fatal(FARGS
,"Error %d initializing PME",status
);
629 if (integrator
[inputrec
->eI
].func
== do_md
)
631 /* Turn on signal handling on all nodes */
633 * (A user signal from the PME nodes (if any)
634 * is communicated to the PP nodes.
636 if (getenv("GMX_NO_TERM") == NULL
)
640 fprintf(debug
,"Installing signal handler for SIGTERM\n");
642 signal(SIGTERM
,signal_handler
);
645 if (getenv("GMX_NO_USR1") == NULL
)
649 fprintf(debug
,"Installing signal handler for SIGUSR1\n");
651 signal(SIGUSR1
,signal_handler
);
656 if (cr
->duty
& DUTY_PP
)
658 if (inputrec
->ePull
!= epullNO
)
660 /* Initialize pull code */
661 init_pull(fplog
,inputrec
,nfile
,fnm
,mtop
,cr
,oenv
,
662 EI_DYNAMICS(inputrec
->eI
) && MASTER(cr
),Flags
);
667 /* Initialize enforced rotation code */
668 init_rot(fplog
,inputrec
,nfile
,fnm
,cr
,box
,state
->x
,oenv
,Flags
);
671 constr
= init_constraints(fplog
,mtop
,inputrec
,ed
,state
,cr
);
673 if (DOMAINDECOMP(cr
))
675 dd_init_bondeds(fplog
,cr
->dd
,mtop
,vsite
,constr
,inputrec
,
676 Flags
& MD_DDBONDCHECK
,fr
->cginfo_mb
);
678 set_dd_parameters(fplog
,cr
->dd
,dlb_scale
,inputrec
,fr
,&ddbox
);
680 setup_dd_grid(fplog
,cr
->dd
);
683 /* Now do whatever the user wants us to do (how flexible...) */
684 integrator
[inputrec
->eI
].func(fplog
,cr
,nfile
,fnm
,
685 oenv
,bVerbose
,bCompact
,
688 nstepout
,inputrec
,mtop
,
690 mdatoms
,nrnb
,wcycle
,ed
,fr
,
691 repl_ex_nst
,repl_ex_seed
,
692 cpt_period
,max_hours
,
696 if (inputrec
->ePull
!= epullNO
)
698 finish_pull(fplog
,inputrec
->pull
);
704 gmx_pmeonly(*pmedata
,cr
,nrnb
,wcycle
,ewaldcoeff
,FALSE
,inputrec
);
707 if (EI_DYNAMICS(inputrec
->eI
) || EI_TPI(inputrec
->eI
))
709 /* Some timing stats */
712 if (runtime
.proc
== 0)
714 runtime
.proc
= runtime
.real
;
723 wallcycle_stop(wcycle
,ewcRUN
);
725 /* Finish up, write some stuff
726 * if rerunMD, don't write last frame again
728 finish_run(fplog
,cr
,ftp2fn(efSTO
,nfile
,fnm
),
729 inputrec
,nrnb
,wcycle
,&runtime
,
730 EI_DYNAMICS(inputrec
->eI
) && !MULTISIM(cr
));
732 /* Does what it says */
733 print_date_and_time(fplog
,cr
->nodeid
,"Finished mdrun",&runtime
);
735 /* Close logfile already here if we were appending to it */
736 if (MASTER(cr
) && (Flags
& MD_APPENDFILES
))
738 gmx_log_close(fplog
);
745 else if(bGotUsr1Signal
)
757 static void md_print_warning(const t_commrec
*cr
,FILE *fplog
,const char *buf
)
761 fprintf(stderr
,"\n%s\n",buf
);
765 fprintf(fplog
,"\n%s\n",buf
);
769 static void check_nst_param(FILE *fplog
,t_commrec
*cr
,
770 const char *desc_nst
,int nst
,
771 const char *desc_p
,int *p
)
775 if (*p
> 0 && *p
% nst
!= 0)
777 /* Round up to the next multiple of nst */
778 *p
= ((*p
)/nst
+ 1)*nst
;
779 sprintf(buf
,"NOTE: %s changes %s to %d\n",desc_nst
,desc_p
,*p
);
780 md_print_warning(cr
,fplog
,buf
);
784 static void reset_all_counters(FILE *fplog
,t_commrec
*cr
,
786 gmx_step_t
*step_rel
,t_inputrec
*ir
,
787 gmx_wallcycle_t wcycle
,t_nrnb
*nrnb
,
788 gmx_runtime_t
*runtime
)
790 char buf
[STRLEN
],sbuf
[22];
792 /* Reset all the counters related to performance over the run */
793 sprintf(buf
,"Step %s: resetting all time and cycle counters\n",
794 gmx_step_str(step
,sbuf
));
795 md_print_warning(cr
,fplog
,buf
);
797 wallcycle_stop(wcycle
,ewcRUN
);
798 wallcycle_reset_all(wcycle
);
799 if (DOMAINDECOMP(cr
))
801 reset_dd_statistics_counters(cr
->dd
);
804 ir
->init_step
+= *step_rel
;
805 ir
->nsteps
-= *step_rel
;
807 wallcycle_start(wcycle
,ewcRUN
);
808 runtime_start(runtime
);
809 print_date_and_time(fplog
,cr
->nodeid
,"Restarted time",runtime
);
812 static int check_nstglobalcomm(FILE *fplog
,t_commrec
*cr
,
813 int nstglobalcomm
,t_inputrec
*ir
)
817 if (!EI_DYNAMICS(ir
->eI
))
822 if (nstglobalcomm
== -1)
824 if (ir
->nstcalcenergy
== 0 && ir
->nstlist
== 0)
827 if (ir
->nstenergy
> 0 && ir
->nstenergy
< nstglobalcomm
)
829 nstglobalcomm
= ir
->nstenergy
;
834 /* We assume that if nstcalcenergy > nstlist,
835 * nstcalcenergy is a multiple of nstlist.
837 if (ir
->nstcalcenergy
== 0 ||
838 (ir
->nstlist
> 0 && ir
->nstlist
< ir
->nstcalcenergy
))
840 nstglobalcomm
= ir
->nstlist
;
844 nstglobalcomm
= ir
->nstcalcenergy
;
850 if (ir
->nstlist
> 0 &&
851 nstglobalcomm
> ir
->nstlist
&& nstglobalcomm
% ir
->nstlist
!= 0)
853 nstglobalcomm
= (nstglobalcomm
/ ir
->nstlist
)*ir
->nstlist
;
854 sprintf(buf
,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm
);
855 md_print_warning(cr
,fplog
,buf
);
857 if (nstglobalcomm
> ir
->nstcalcenergy
)
859 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
860 "nstcalcenergy",&ir
->nstcalcenergy
);
863 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
864 "nstenergy",&ir
->nstenergy
);
866 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
867 "nstlog",&ir
->nstlog
);
870 if (ir
->comm_mode
!= ecmNO
&& ir
->nstcomm
< nstglobalcomm
)
872 sprintf(buf
,"WARNING: Changing nstcomm from %d to %d\n",
873 ir
->nstcomm
,nstglobalcomm
);
874 md_print_warning(cr
,fplog
,buf
);
875 ir
->nstcomm
= nstglobalcomm
;
878 return nstglobalcomm
;
881 static void check_ir_old_tpx_versions(t_commrec
*cr
,FILE *fplog
,
882 t_inputrec
*ir
,gmx_mtop_t
*mtop
)
884 /* Check required for old tpx files */
885 if (IR_TWINRANGE(*ir
) && ir
->nstlist
> 1 &&
886 ir
->nstcalcenergy
% ir
->nstlist
!= 0)
888 md_print_warning(cr
,fplog
,"Old tpr file with twin-range settings: modifiying energy calculation and/or T/P-coupling frequencies");
890 if (gmx_mtop_ftype_count(mtop
,F_CONSTR
) +
891 gmx_mtop_ftype_count(mtop
,F_CONSTRNC
) > 0 &&
892 ir
->eConstrAlg
== econtSHAKE
)
894 md_print_warning(cr
,fplog
,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect");
895 if (ir
->epc
!= epcNO
)
897 gmx_fatal(FARGS
,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
900 check_nst_param(fplog
,cr
,"nstlist",ir
->nstlist
,
901 "nstcalcenergy",&ir
->nstcalcenergy
);
903 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
904 "nstenergy",&ir
->nstenergy
);
905 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
906 "nstlog",&ir
->nstlog
);
907 if (ir
->efep
!= efepNO
)
909 check_nst_param(fplog
,cr
,"nstdhdl",ir
->nstdhdl
,
910 "nstenergy",&ir
->nstenergy
);
915 bool bGStatEveryStep
;
917 gmx_step_t step_nscheck
;
928 static void reset_nlistheuristics(gmx_nlheur_t
*nlh
,gmx_step_t step
)
932 nlh
->step_nscheck
= step
;
935 static void init_nlistheuristics(gmx_nlheur_t
*nlh
,
936 bool bGStatEveryStep
,gmx_step_t step
)
938 nlh
->bGStatEveryStep
= bGStatEveryStep
;
945 reset_nlistheuristics(nlh
,step
);
948 static void update_nliststatistics(gmx_nlheur_t
*nlh
,gmx_step_t step
)
951 char sbuf
[22],sbuf2
[22];
953 /* Determine the neighbor list life time */
954 nl_lt
= step
- nlh
->step_ns
;
957 fprintf(debug
,"%d atoms beyond ns buffer, updating neighbor list after %s steps\n",nlh
->nabnsb
,gmx_step_str(nl_lt
,sbuf
));
961 nlh
->s2
+= nl_lt
*nl_lt
;
962 nlh
->ab
+= nlh
->nabnsb
;
963 if (nlh
->lt_runav
== 0)
965 nlh
->lt_runav
= nl_lt
;
966 /* Initialize the fluctuation average
967 * such that at startup we check after 0 steps.
969 nlh
->lt_runav2
= sqr(nl_lt
/2.0);
971 /* Running average with 0.9 gives an exp. history of 9.5 */
972 nlh
->lt_runav2
= 0.9*nlh
->lt_runav2
+ 0.1*sqr(nlh
->lt_runav
- nl_lt
);
973 nlh
->lt_runav
= 0.9*nlh
->lt_runav
+ 0.1*nl_lt
;
974 if (nlh
->bGStatEveryStep
)
976 /* Always check the nlist validity */
977 nlh
->step_nscheck
= step
;
981 /* We check after: <life time> - 2*sigma
982 * The factor 2 is quite conservative,
983 * but we assume that with nstlist=-1 the user
984 * prefers exact integration over performance.
986 nlh
->step_nscheck
= step
987 + (int)(nlh
->lt_runav
- 2.0*sqrt(nlh
->lt_runav2
)) - 1;
991 fprintf(debug
,"nlist life time %s run av. %4.1f sig %3.1f check %s check with -gcom %d\n",
992 gmx_step_str(nl_lt
,sbuf
),nlh
->lt_runav
,sqrt(nlh
->lt_runav2
),
993 gmx_step_str(nlh
->step_nscheck
-step
+1,sbuf2
),
994 (int)(nlh
->lt_runav
- 2.0*sqrt(nlh
->lt_runav2
)));
998 static void set_nlistheuristics(gmx_nlheur_t
*nlh
,bool bReset
,gmx_step_t step
)
1004 reset_nlistheuristics(nlh
,step
);
1008 update_nliststatistics(nlh
,step
);
1011 nlh
->step_ns
= step
;
1012 /* Initialize the cumulative coordinate scaling matrix */
1013 clear_mat(nlh
->scale_tot
);
1014 for(d
=0; d
<DIM
; d
++)
1016 nlh
->scale_tot
[d
][d
] = 1.0;
1020 double do_md(FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
1021 const output_env_t oenv
, bool bVerbose
,bool bCompact
,
1023 gmx_vsite_t
*vsite
,gmx_constr_t constr
,
1024 int stepout
,t_inputrec
*ir
,
1025 gmx_mtop_t
*top_global
,
1027 t_state
*state_global
,
1029 t_nrnb
*nrnb
,gmx_wallcycle_t wcycle
,
1030 gmx_edsam_t ed
,t_forcerec
*fr
,
1031 int repl_ex_nst
,int repl_ex_seed
,
1032 real cpt_period
,real max_hours
,
1033 unsigned long Flags
,
1034 gmx_runtime_t
*runtime
)
1036 int fp_trn
=0,fp_xtc
=0;
1037 ener_file_t fp_ene
=NULL
;
1038 gmx_step_t step
,step_rel
;
1040 FILE *fp_dhdl
=NULL
,*fp_field
=NULL
;
1043 bool bGStatEveryStep
,bGStat
,bNstEner
,bCalcEner
;
1044 bool bNS
,bNStList
,bSimAnn
,bStopCM
,bRerunMD
,bNotLastFrame
=FALSE
,
1045 bFirstStep
,bStateFromTPX
,bLastStep
,bBornRadii
;
1047 bool bNEMD
,do_ene
,do_log
,do_verbose
,bRerunWarnNoV
=TRUE
,
1048 bForceUpdate
=FALSE
,bX
,bV
,bF
,bXTC
,bCPT
;
1051 tensor force_vir
,shake_vir
,total_vir
,pres
,ekin
;
1056 t_trxframe rerun_fr
;
1057 gmx_repl_ex_t repl_ex
=NULL
;
1059 /* Booleans (disguised as a reals) to checkpoint and terminate mdrun */
1060 real chkpt
=0,terminate
=0,terminate_now
=0;
1062 gmx_localtop_t
*top
;
1063 t_mdebin
*mdebin
=NULL
;
1064 t_state
*state
=NULL
;
1065 rvec
*f_global
=NULL
;
1068 gmx_enerdata_t
*enerd
;
1070 gmx_global_stat_t gstat
;
1071 gmx_update_t upd
=NULL
;
1072 t_graph
*graph
=NULL
;
1075 gmx_groups_t
*groups
;
1076 gmx_ekindata_t
*ekind
;
1077 gmx_shellfc_t shellfc
;
1078 int count
,nconverged
=0;
1082 bool bTCR
=FALSE
,bConverged
=TRUE
,bOK
,bSumEkinhOld
,bExchanged
;
1084 real temp0
,mu_aver
=0,dvdl
;
1086 atom_id
*grpindex
=NULL
;
1088 t_coupl_rec
*tcr
=NULL
;
1089 rvec
*xcopy
=NULL
,*vcopy
=NULL
;
1090 matrix boxcopy
,lastbox
;
1092 int reset_counters
=-1;
1093 char sbuf
[22],sbuf2
[22];
1094 bool bHandledSignal
=FALSE
;
1096 /* Temporary addition for FAHCORE checkpointing */
1101 /* Check for special mdrun options */
1102 bRerunMD
= (Flags
& MD_RERUN
);
1103 bIonize
= (Flags
& MD_IONIZE
);
1104 bFFscan
= (Flags
& MD_FFSCAN
);
1105 bAppend
= (Flags
& MD_APPENDFILES
);
1109 /* Since we don't know if the frames read are related in any way,
1110 * rebuild the neighborlist at every step.
1113 ir
->nstcalcenergy
= 1;
1117 check_ir_old_tpx_versions(cr
,fplog
,ir
,top_global
);
1119 nstglobalcomm
= check_nstglobalcomm(fplog
,cr
,nstglobalcomm
,ir
);
1120 bGStatEveryStep
= (nstglobalcomm
== 1);
1122 if (!bGStatEveryStep
&& ir
->nstlist
== -1 && fplog
!= NULL
)
1125 "To reduce the energy communication with nstlist = -1\n"
1126 "the neighbor list validity should not be checked at every step,\n"
1127 "this means that exact integration is not guaranteed.\n"
1128 "The neighbor list validity is checked after:\n"
1129 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
1130 "In most cases this will result in exact integration.\n"
1131 "This reduces the energy communication by a factor of 2 to 3.\n"
1132 "If you want less energy communication, set nstlist > 3.\n\n");
1135 if (bRerunMD
|| bFFscan
)
1139 groups
= &top_global
->groups
;
1141 /* Initial values */
1142 init_md(fplog
,cr
,ir
,oenv
,&t
,&t0
,&state_global
->lambda
,&lam0
,
1143 nrnb
,top_global
,&upd
,
1144 nfile
,fnm
,&fp_trn
,&fp_xtc
,&fp_ene
,&fn_cpt
,
1145 &fp_dhdl
,&fp_field
,&mdebin
,
1146 force_vir
,shake_vir
,mu_tot
,&bNEMD
,&bSimAnn
,&vcm
,Flags
);
1148 /* Energy terms and groups */
1150 init_enerdata(top_global
->groups
.grps
[egcENER
].nr
,ir
->n_flambda
,enerd
);
1151 if (DOMAINDECOMP(cr
))
1157 snew(f
,top_global
->natoms
);
1160 /* Kinetic energy data */
1162 init_ekindata(fplog
,top_global
,&(ir
->opts
),ekind
);
1163 /* Copy the cos acceleration to the groups struct */
1164 ekind
->cosacc
.cos_accel
= ir
->cos_accel
;
1166 gstat
= global_stat_init(ir
);
1169 /* Check for polarizable models and flexible constraints */
1170 shellfc
= init_shell_flexcon(fplog
,
1171 top_global
,n_flexible_constraints(constr
),
1172 (ir
->bContinuation
||
1173 (DOMAINDECOMP(cr
) && !MASTER(cr
))) ?
1174 NULL
: state_global
->x
);
1179 tMPI_Thread_mutex_lock(&box_mutex
);
1181 set_deform_reference_box(upd
,init_step_tpx
,box_tpx
);
1183 tMPI_Thread_mutex_unlock(&box_mutex
);
1188 double io
= compute_io(ir
,top_global
->natoms
,groups
,mdebin
->ebin
->nener
,1);
1189 if ((io
> 2000) && MASTER(cr
))
1191 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
1195 if (DOMAINDECOMP(cr
)) {
1196 top
= dd_init_local_top(top_global
);
1199 dd_init_local_state(cr
->dd
,state_global
,state
);
1201 if (DDMASTER(cr
->dd
) && ir
->nstfout
) {
1202 snew(f_global
,state_global
->natoms
);
1206 /* Initialize the particle decomposition and split the topology */
1207 top
= split_system(fplog
,top_global
,ir
,cr
);
1209 pd_cg_range(cr
,&fr
->cg0
,&fr
->hcg
);
1210 pd_at_range(cr
,&a0
,&a1
);
1212 top
= gmx_mtop_generate_local_top(top_global
,ir
);
1215 a1
= top_global
->natoms
;
1218 state
= partdec_init_local_state(cr
,state_global
);
1221 atoms2md(top_global
,ir
,0,NULL
,a0
,a1
-a0
,mdatoms
);
1224 set_vsite_top(vsite
,top
,mdatoms
,cr
);
1227 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
) {
1228 graph
= mk_graph(fplog
,&(top
->idef
),0,top_global
->natoms
,FALSE
,FALSE
);
1232 make_local_shells(cr
,mdatoms
,shellfc
);
1235 if (ir
->pull
&& PAR(cr
)) {
1236 dd_make_local_pull_groups(NULL
,ir
->pull
,mdatoms
);
1240 if (DOMAINDECOMP(cr
))
1242 /* Distribute the charge groups over the nodes from the master node */
1243 dd_partition_system(fplog
,ir
->init_step
,cr
,TRUE
,1,
1244 state_global
,top_global
,ir
,
1245 state
,&f
,mdatoms
,top
,fr
,
1246 vsite
,shellfc
,constr
,
1250 /* If not DD, copy gb data */
1251 if(ir
->implicit_solvent
&& !DOMAINDECOMP(cr
))
1253 make_local_gb(cr
,fr
->born
,ir
->gb_algorithm
);
1256 update_mdatoms(mdatoms
,state
->lambda
);
1260 /* Update mdebin with energy history if appending to output files */
1261 if ( Flags
& MD_APPENDFILES
)
1263 restore_energyhistory_from_state(mdebin
,&state_global
->enerhist
);
1265 /* Set the initial energy history in state to zero by updating once */
1266 update_energyhistory(&state_global
->enerhist
,mdebin
);
1270 if ((state
->flags
& (1<<estLD_RNG
)) && (Flags
& MD_READ_RNG
)) {
1271 /* Set the random state if we read a checkpoint file */
1272 set_stochd_state(upd
,state
);
1275 /* Initialize constraints */
1277 if (!DOMAINDECOMP(cr
))
1278 set_constraints(constr
,top
,ir
,mdatoms
,cr
);
1281 /* Check whether we have to GCT stuff */
1282 bTCR
= ftp2bSet(efGCT
,nfile
,fnm
);
1285 fprintf(stderr
,"Will do General Coupling Theory!\n");
1287 gnx
= top_global
->mols
.nr
;
1289 for(i
=0; (i
<gnx
); i
++) {
1294 if (repl_ex_nst
> 0 && MASTER(cr
))
1295 repl_ex
= init_replica_exchange(fplog
,cr
->ms
,state_global
,ir
,
1296 repl_ex_nst
,repl_ex_seed
);
1299 if (!ir
->bContinuation
&& !bRerunMD
)
1301 if (mdatoms
->cFREEZE
&& (state
->flags
& (1<<estV
)))
1303 /* Set the velocities of frozen particles to zero */
1304 for(i
=mdatoms
->start
; i
<mdatoms
->start
+mdatoms
->homenr
; i
++)
1306 for(m
=0; m
<DIM
; m
++)
1308 if (ir
->opts
.nFreeze
[mdatoms
->cFREEZE
[i
]][m
])
1318 /* Constrain the initial coordinates and velocities */
1319 do_constrain_first(fplog
,constr
,ir
,mdatoms
,state
,
1320 graph
,cr
,nrnb
,fr
,&top
->idef
);
1324 /* Construct the virtual sites for the initial configuration */
1325 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,NULL
,
1326 top
->idef
.iparams
,top
->idef
.il
,
1327 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
1333 if (Flags
& MD_READ_EKIN
)
1335 restore_ekinstate_from_state(cr
,ekind
,&state_global
->ekinstate
);
1339 /* Compute initial EKin for all.. */
1340 calc_ke_part(state
,&(ir
->opts
),mdatoms
,ekind
,nrnb
);
1345 GMX_MPE_LOG(ev_global_stat_start
);
1347 global_stat(fplog
,gstat
,cr
,enerd
,force_vir
,shake_vir
,mu_tot
,
1348 ir
,ekind
,FALSE
,constr
,vcm
,NULL
,NULL
,&terminate
,
1351 GMX_MPE_LOG(ev_global_stat_finish
);
1356 /* Calculate the initial half step temperature */
1357 temp0
= sum_ekin(TRUE
,&(ir
->opts
),ekind
,ekin
,NULL
);
1361 /* Initiate data for the special cases */
1363 snew(xcopy
,state
->natoms
);
1364 snew(vcopy
,state
->natoms
);
1365 for(ii
=0; (ii
<state
->natoms
); ii
++) {
1366 copy_rvec(state
->x
[ii
],xcopy
[ii
]);
1367 copy_rvec(state
->v
[ii
],vcopy
[ii
]);
1369 copy_mat(state
->box
,boxcopy
);
1374 if (constr
&& !ir
->bContinuation
&& ir
->eConstrAlg
== econtLINCS
)
1377 "RMS relative constraint deviation after constraining: %.2e\n",
1378 constr_rmsd(constr
,FALSE
));
1380 fprintf(fplog
,"Initial temperature: %g K\n",temp0
);
1383 fprintf(stderr
,"starting md rerun '%s', reading coordinates from"
1384 " input trajectory '%s'\n\n",
1385 *(top_global
->name
),opt2fn("-rerun",nfile
,fnm
));
1388 fprintf(stderr
,"Calculated time to finish depends on nsteps from "
1389 "run input file,\nwhich may not correspond to the time "
1390 "needed to process input trajectory.\n\n");
1395 fprintf(stderr
,"starting mdrun '%s'\n",
1396 *(top_global
->name
));
1397 if (ir
->init_step
> 0)
1399 fprintf(stderr
,"%s steps, %8.1f ps (continuing from step %s, %8.1f ps).\n",
1400 gmx_step_str(ir
->init_step
+ir
->nsteps
,sbuf
),
1401 (ir
->init_step
+ir
->nsteps
)*ir
->delta_t
,
1402 gmx_step_str(ir
->init_step
,sbuf2
),
1403 ir
->init_step
*ir
->delta_t
);
1407 fprintf(stderr
,"%s steps, %8.1f ps.\n",
1408 gmx_step_str(ir
->nsteps
,sbuf
),ir
->nsteps
*ir
->delta_t
);
1411 fprintf(fplog
,"\n");
1414 /* Set and write start time */
1415 runtime_start(runtime
);
1416 print_date_and_time(fplog
,cr
->nodeid
,"Started mdrun",runtime
);
1417 wallcycle_start(wcycle
,ewcRUN
);
1419 fprintf(fplog
,"\n");
1421 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
1423 chkpt_ret
=fcCheckPointParallel( cr
->nodeid
,
1425 if ( chkpt_ret
== 0 )
1426 gmx_fatal( 3,__FILE__
,__LINE__
, "Checkpoint error on step %d\n", 0 );
1430 /***********************************************************
1432 * Loop over MD steps
1434 ************************************************************/
1436 /* if rerunMD then read coordinates and velocities from input trajectory */
1439 if (getenv("GMX_FORCE_UPDATE"))
1441 bForceUpdate
= TRUE
;
1444 bNotLastFrame
= read_first_frame(oenv
,&status
,
1445 opt2fn("-rerun",nfile
,fnm
),
1446 &rerun_fr
,TRX_NEED_X
| TRX_READ_V
);
1447 if (rerun_fr
.natoms
!= top_global
->natoms
)
1450 "Number of atoms in trajectory (%d) does not match the "
1451 "run input file (%d)\n",
1452 rerun_fr
.natoms
,top_global
->natoms
);
1454 if (ir
->ePBC
!= epbcNONE
)
1458 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
);
1460 if (max_cutoff2(ir
->ePBC
,rerun_fr
.box
) < sqr(fr
->rlistlong
))
1462 gmx_fatal(FARGS
,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr
.step
,rerun_fr
.time
);
1465 /* Set the shift vectors.
1466 * Necessary here when have a static box different from the tpr box.
1468 calc_shifts(rerun_fr
.box
,fr
->shift_vec
);
1472 /* loop over MD steps or if rerunMD to end of input trajectory */
1474 /* Skip the first Nose-Hoover integration when we get the state from tpx */
1475 bStateFromTPX
= !opt2bSet("-cpi",nfile
,fnm
);
1477 bSumEkinhOld
= FALSE
;
1480 step
= ir
->init_step
;
1483 if (ir
->nstlist
== -1)
1485 init_nlistheuristics(&nlh
,bGStatEveryStep
,step
);
1488 bLastStep
= (bRerunMD
|| step_rel
> ir
->nsteps
);
1489 while (!bLastStep
|| (bRerunMD
&& bNotLastFrame
)) {
1491 wallcycle_start(wcycle
,ewcSTEP
);
1493 GMX_MPE_LOG(ev_timestep1
);
1496 if (rerun_fr
.bStep
) {
1497 step
= rerun_fr
.step
;
1498 step_rel
= step
- ir
->init_step
;
1505 bLastStep
= (step_rel
== ir
->nsteps
);
1507 t
= t0
+ step
*ir
->delta_t
;
1510 if (ir
->efep
!= efepNO
)
1512 if (bRerunMD
&& rerun_fr
.bLambda
&& (ir
->delta_lambda
!=0))
1514 state_global
->lambda
= rerun_fr
.lambda
;
1518 state_global
->lambda
= lam0
+ step
*ir
->delta_lambda
;
1520 state
->lambda
= state_global
->lambda
;
1521 bDoDHDL
= do_per_step(step
,ir
->nstdhdl
);
1526 update_annealing_target_temp(&(ir
->opts
),t
);
1531 if (!(DOMAINDECOMP(cr
) && !MASTER(cr
)))
1533 for(i
=0; i
<state_global
->natoms
; i
++)
1535 copy_rvec(rerun_fr
.x
[i
],state_global
->x
[i
]);
1539 for(i
=0; i
<state_global
->natoms
; i
++)
1541 copy_rvec(rerun_fr
.v
[i
],state_global
->v
[i
]);
1546 for(i
=0; i
<state_global
->natoms
; i
++)
1548 clear_rvec(state_global
->v
[i
]);
1552 fprintf(stderr
,"\nWARNING: Some frames do not contain velocities.\n"
1553 " Ekin, temperature and pressure are incorrect,\n"
1554 " the virial will be incorrect when constraints are present.\n"
1556 bRerunWarnNoV
= FALSE
;
1560 copy_mat(rerun_fr
.box
,state_global
->box
);
1561 copy_mat(state_global
->box
,state
->box
);
1563 if (vsite
&& (Flags
& MD_RERUN_VSITE
))
1565 if (DOMAINDECOMP(cr
))
1567 gmx_fatal(FARGS
,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
1571 /* Following is necessary because the graph may get out of sync
1572 * with the coordinates if we only have every N'th coordinate set
1574 mk_mshift(fplog
,graph
,fr
->ePBC
,state
->box
,state
->x
);
1575 shift_self(graph
,state
->box
,state
->x
);
1577 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,state
->v
,
1578 top
->idef
.iparams
,top
->idef
.il
,
1579 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
1582 unshift_self(graph
,state
->box
,state
->x
);
1587 /* Stop Center of Mass motion */
1588 bStopCM
= (ir
->comm_mode
!= ecmNO
&& do_per_step(step
,ir
->nstcomm
));
1590 /* Copy back starting coordinates in case we're doing a forcefield scan */
1593 for(ii
=0; (ii
<state
->natoms
); ii
++)
1595 copy_rvec(xcopy
[ii
],state
->x
[ii
]);
1596 copy_rvec(vcopy
[ii
],state
->v
[ii
]);
1598 copy_mat(boxcopy
,state
->box
);
1603 /* for rerun MD always do Neighbour Searching */
1604 bNS
= (bFirstStep
|| ir
->nstlist
!= 0);
1609 /* Determine whether or not to do Neighbour Searching and LR */
1610 bNStList
= (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0);
1612 bNS
= (bFirstStep
|| bExchanged
|| bNStList
||
1613 (ir
->nstlist
== -1 && nlh
.nabnsb
> 0));
1615 if (bNS
&& ir
->nstlist
== -1)
1617 set_nlistheuristics(&nlh
,bFirstStep
|| bExchanged
,step
);
1621 if (terminate_now
> 0 || (terminate_now
< 0 && bNS
))
1626 /* Determine whether or not to update the Born radii if doing GB */
1627 bBornRadii
=bFirstStep
;
1628 if(ir
->implicit_solvent
&& (step
% ir
->nstgbradii
==0))
1633 do_log
= do_per_step(step
,ir
->nstlog
) || bFirstStep
|| bLastStep
;
1634 do_verbose
= bVerbose
&&
1635 (step
% stepout
== 0 || bFirstStep
|| bLastStep
);
1637 if (bNS
&& !(bFirstStep
&& ir
->bContinuation
&& !bRerunMD
))
1641 bMasterState
= TRUE
;
1645 bMasterState
= FALSE
;
1646 /* Correct the new box if it is too skewed */
1647 if (DYNAMIC_BOX(*ir
))
1649 if (correct_box(fplog
,step
,state
->box
,graph
))
1651 bMasterState
= TRUE
;
1654 if (DOMAINDECOMP(cr
) && bMasterState
)
1656 dd_collect_state(cr
->dd
,state
,state_global
);
1660 if (DOMAINDECOMP(cr
))
1662 /* Repartition the domain decomposition */
1663 wallcycle_start(wcycle
,ewcDOMDEC
);
1664 dd_partition_system(fplog
,step
,cr
,
1665 bMasterState
,nstglobalcomm
,
1666 state_global
,top_global
,ir
,
1667 state
,&f
,mdatoms
,top
,fr
,
1668 vsite
,shellfc
,constr
,
1669 nrnb
,wcycle
,do_verbose
);
1670 wallcycle_stop(wcycle
,ewcDOMDEC
);
1674 if (MASTER(cr
) && do_log
&& !bFFscan
)
1676 print_ebin_header(fplog
,step
,t
,state
->lambda
);
1679 if (ir
->efep
!= efepNO
)
1681 update_mdatoms(mdatoms
,state
->lambda
);
1684 if (bRerunMD
&& rerun_fr
.bV
)
1686 /* We need the kinetic energy at minus the half step for determining
1687 * the full step kinetic energy and possibly for T-coupling.
1689 calc_ke_part(state
,&(ir
->opts
),mdatoms
,ekind
,nrnb
);
1692 global_stat(fplog
,gstat
,cr
,enerd
,force_vir
,shake_vir
,mu_tot
,
1693 ir
,ekind
,FALSE
,constr
,vcm
,NULL
,NULL
,&terminate
,
1696 sum_ekin(FALSE
,&(ir
->opts
),ekind
,ekin
,NULL
);
1698 clear_mat(force_vir
);
1700 /* Ionize the atoms if necessary */
1703 ionize(fplog
,oenv
,mdatoms
,top_global
,t
,ir
,state
->x
,state
->v
,
1704 mdatoms
->start
,mdatoms
->start
+mdatoms
->homenr
,state
->box
,cr
);
1707 /* Update force field in ffscan program */
1710 if (update_forcefield(fplog
,
1712 mdatoms
->nr
,state
->x
,state
->box
)) {
1713 if (gmx_parallel_env())
1721 GMX_MPE_LOG(ev_timestep2
);
1723 if ((bNS
|| bLastStep
) && (step
> ir
->init_step
) && !bRerunMD
)
1725 bCPT
= (chkpt
> 0 || (bLastStep
&& (Flags
& MD_CONFOUT
)));
1736 /* Determine the pressure:
1737 * always when we want exact averages in the energy file,
1738 * at ns steps when we have pressure coupling,
1739 * otherwise only at energy output steps (set below).
1741 bNstEner
= (bGStatEveryStep
|| do_per_step(step
,ir
->nstcalcenergy
));
1742 bCalcEner
= bNstEner
;
1744 /* Do we need global communication ? */
1745 bGStat
= (bCalcEner
|| bStopCM
||
1746 (ir
->nstlist
== -1 && !bRerunMD
&& step
>= nlh
.step_nscheck
));
1748 do_ene
= (do_per_step(step
,ir
->nstenergy
) || bLastStep
);
1750 if (do_ene
|| do_log
)
1756 force_flags
= (GMX_FORCE_STATECHANGED
|
1757 GMX_FORCE_ALLFORCES
|
1758 (bNStList
? GMX_FORCE_DOLR
: 0) |
1760 (bCalcEner
? GMX_FORCE_VIRIAL
: 0) |
1761 (bDoDHDL
? GMX_FORCE_DHDL
: 0));
1765 /* Now is the time to relax the shells */
1766 count
=relax_shell_flexcon(fplog
,cr
,bVerbose
,bFFscan
? step
+1 : step
,
1768 bStopCM
,top
,top_global
,
1770 state
,f
,force_vir
,mdatoms
,
1771 nrnb
,wcycle
,graph
,groups
,
1772 shellfc
,fr
,bBornRadii
,t
,mu_tot
,
1773 state
->natoms
,&bConverged
,vsite
,
1784 /* The coordinates (x) are shifted (to get whole molecules)
1786 * This is parallellized as well, and does communication too.
1787 * Check comments in sim_util.c
1790 do_force(fplog
,cr
,ir
,step
,nrnb
,wcycle
,top
,top_global
,groups
,
1791 state
->box
,state
->x
,&state
->hist
,
1792 f
,force_vir
,mdatoms
,enerd
,fcd
,
1793 state
->lambda
,graph
,
1794 fr
,vsite
,mu_tot
,t
,fp_field
,ed
,bBornRadii
,
1795 (bNS
? GMX_FORCE_NS
: 0) | force_flags
);
1798 GMX_BARRIER(cr
->mpi_comm_mygroup
);
1802 mu_aver
= calc_mu_aver(cr
,state
->x
,mdatoms
->chargeA
,
1803 mu_tot
,&top_global
->mols
,mdatoms
,gnx
,grpindex
);
1806 if (bTCR
&& bFirstStep
)
1808 tcr
=init_coupling(fplog
,nfile
,fnm
,cr
,fr
,mdatoms
,&(top
->idef
));
1809 fprintf(fplog
,"Done init_coupling\n");
1813 /* Now we have the energies and forces corresponding to the
1814 * coordinates at time t. We must output all of this before
1816 * for RerunMD t is read from input trajectory
1818 GMX_MPE_LOG(ev_output_start
);
1820 bX
= do_per_step(step
,ir
->nstxout
);
1821 bV
= do_per_step(step
,ir
->nstvout
);
1822 bF
= do_per_step(step
,ir
->nstfout
);
1823 bXTC
= do_per_step(step
,ir
->nstxtcout
);
1827 fcReportProgress( ir
->nsteps
, step
);
1829 bX
= bX
|| bLastStep
; /*enforce writing positions and velocities
1831 bV
= bV
|| bLastStep
;
1833 int nthreads
=(cr
->nthreads
==0 ? 1 : cr
->nthreads
);
1834 int nnodes
=(cr
->nnodes
==0 ? 1 : cr
->nnodes
);
1837 /*Gromacs drives checkpointing; no ||
1838 fcCheckPointPendingThreads(cr->nodeid,
1840 /* sync bCPT and fc record-keeping */
1841 if (bCPT
&& MASTER(cr
))
1842 fcRequestCheckPoint();
1847 if (bX
|| bV
|| bF
|| bXTC
|| bCPT
)
1849 wallcycle_start(wcycle
,ewcTRAJ
);
1852 if (state
->flags
& (1<<estLD_RNG
))
1854 get_stochd_state(upd
,state
);
1860 state_global
->ekinstate
.bUpToDate
= FALSE
;
1864 update_ekinstate(&state_global
->ekinstate
,ekind
);
1865 state_global
->ekinstate
.bUpToDate
= TRUE
;
1867 update_energyhistory(&state_global
->enerhist
,mdebin
);
1870 write_traj(fplog
,cr
,fp_trn
,bX
,bV
,bF
,fp_xtc
,bXTC
,ir
->xtcprec
,
1871 fn_cpt
,bCPT
,top_global
,ir
->eI
,ir
->simulation_part
,
1872 step
,t
,state
,state_global
,f
,f_global
,&n_xtc
,&x_xtc
);
1874 if (bLastStep
&& step_rel
== ir
->nsteps
&&
1875 (Flags
& MD_CONFOUT
) && MASTER(cr
) &&
1876 !bRerunMD
&& !bFFscan
)
1878 /* x and v have been collected in write_traj */
1879 fprintf(stderr
,"\nWriting final coordinates.\n");
1880 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
&&
1883 /* Make molecules whole only for confout writing */
1884 do_pbc_mtop(fplog
,ir
->ePBC
,state
->box
,top_global
,state_global
->x
);
1886 write_sto_conf_mtop(ftp2fn(efSTO
,nfile
,fnm
),
1887 *top_global
->name
,top_global
,
1888 state_global
->x
,state_global
->v
,
1889 ir
->ePBC
,state
->box
);
1892 wallcycle_stop(wcycle
,ewcTRAJ
);
1894 GMX_MPE_LOG(ev_output_finish
);
1896 clear_mat(shake_vir
);
1898 /* Box is changed in update() when we do pressure coupling,
1899 * but we should still use the old box for energy corrections and when
1900 * writing it to the energy file, so it matches the trajectory files for
1901 * the same timestep above. Make a copy in a separate array.
1903 copy_mat(state
->box
,lastbox
);
1906 GMX_MPE_LOG(ev_update_start
);
1907 /* This is also parallellized, but check code in update.c */
1908 /* bOK = update(nsb->natoms,START(nsb),HOMENR(nsb),step,state->lambda,&ener[F_DVDL], */
1910 if (!bRerunMD
|| rerun_fr
.bV
|| bForceUpdate
)
1912 wallcycle_start(wcycle
,ewcUPDATE
);
1914 /* We can only do Berendsen coupling after we have summed
1915 * the kinetic energy or virial.
1916 * Since the happens in global_state after update,
1917 * we should only do it at step % nstlist = 1
1918 * with bGStatEveryStep=FALSE.
1920 update(fplog
,step
,&dvdl
,ir
,mdatoms
,state
,graph
,
1921 f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
1922 &top
->idef
,ekind
,ir
->nstlist
==-1 ? &nlh
.scale_tot
: NULL
,
1923 cr
,nrnb
,wcycle
,upd
,constr
,bCalcEner
,shake_vir
,
1924 bNEMD
,bFirstStep
&& bStateFromTPX
);
1925 if (fr
->bSepDVDL
&& fplog
&& do_log
)
1927 fprintf(fplog
,sepdvdlformat
,"Constraint",0.0,dvdl
);
1929 enerd
->term
[F_DHDL_CON
] += dvdl
;
1930 wallcycle_stop(wcycle
,ewcUPDATE
);
1934 /* Need to unshift here */
1935 unshift_self(graph
,state
->box
,state
->x
);
1938 GMX_BARRIER(cr
->mpi_comm_mygroup
);
1939 GMX_MPE_LOG(ev_update_finish
);
1941 if (!bOK
&& !bFFscan
)
1943 gmx_fatal(FARGS
,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1948 wallcycle_start(wcycle
,ewcVSITECONSTR
);
1951 shift_self(graph
,state
->box
,state
->x
);
1953 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,state
->v
,
1954 top
->idef
.iparams
,top
->idef
.il
,
1955 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
1959 unshift_self(graph
,state
->box
,state
->x
);
1961 wallcycle_stop(wcycle
,ewcVSITECONSTR
);
1964 /* Non-equilibrium MD: this is parallellized,
1965 * but only does communication when there really is NEMD.
1967 if (PAR(cr
) && bNEMD
)
1969 accumulate_u(cr
,&(ir
->opts
),ekind
);
1973 calc_ke_part(state
,&(ir
->opts
),mdatoms
,ekind
,nrnb
);
1975 /* since we use the new coordinates in calc_ke_part_visc, we should use
1976 * the new box too. Still, won't this be offset by one timestep in the
1977 * energy file? / EL 20040121
1981 /* Calculate center of mass velocity if necessary, also parallellized */
1982 if (bStopCM
&& !bFFscan
&& !bRerunMD
)
1984 calc_vcm_grp(fplog
,mdatoms
->start
,mdatoms
->homenr
,mdatoms
,
1985 state
->x
,state
->v
,vcm
);
1988 /* Determine the wallclock run time up till now */
1989 run_time
= (double)time(NULL
) - (double)runtime
->real
;
1991 /* Check whether everything is still allright */
1992 if ((bGotTermSignal
|| bGotUsr1Signal
) && !bHandledSignal
)
1994 if (bGotTermSignal
|| ir
->nstlist
== 0)
2004 terminate_now
= terminate
;
2009 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2010 bGotTermSignal
? "TERM" : "USR1",
2011 terminate
==-1 ? "NS " : "");
2015 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2016 bGotTermSignal
? "TERM" : "USR1",
2017 terminate
==-1 ? "NS " : "");
2019 bHandledSignal
=TRUE
;
2021 else if (MASTER(cr
) && (bNS
|| ir
->nstlist
<= 0) &&
2022 (max_hours
> 0 && run_time
> max_hours
*60.0*60.0*0.99) &&
2025 /* Signal to terminate the run */
2026 terminate
= (ir
->nstlist
== 0 ? 1 : -1);
2029 terminate_now
= terminate
;
2033 fprintf(fplog
,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step
,sbuf
),max_hours
*0.99);
2035 fprintf(stderr
, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step
,sbuf
),max_hours
*0.99);
2038 if (ir
->nstlist
== -1 && !bRerunMD
)
2040 /* When bGStatEveryStep=FALSE, global_stat is only called
2041 * when we check the atom displacements, not at NS steps.
2042 * This means that also the bonded interaction count check is not
2043 * performed immediately after NS. Therefore a few MD steps could
2044 * be performed with missing interactions.
2045 * But wrong energies are never written to file,
2046 * since energies are only written after global_stat
2049 if (step
>= nlh
.step_nscheck
)
2051 nlh
.nabnsb
= natoms_beyond_ns_buffer(ir
,fr
,&top
->cgs
,
2052 nlh
.scale_tot
,state
->x
);
2056 /* This is not necessarily true,
2057 * but step_nscheck is determined quite conservatively.
2063 /* In parallel we only have to check for checkpointing in steps
2064 * where we do global communication,
2065 * otherwise the other nodes don't know.
2067 if (MASTER(cr
) && ((bGStat
|| !PAR(cr
)) &&
2070 run_time
>= nchkpt
*cpt_period
*60.0)))
2081 /* We will not sum ekinh_old,
2082 * so signal that we still have to do it.
2084 bSumEkinhOld
= TRUE
;
2090 wallcycle_start(wcycle
,ewcMoveE
);
2091 /* Globally (over all NODEs) sum energy, virial etc.
2092 * This includes communication
2094 global_stat(fplog
,gstat
,cr
,enerd
,force_vir
,shake_vir
,mu_tot
,
2095 ir
,ekind
,bSumEkinhOld
,constr
,vcm
,
2096 ir
->nstlist
==-1 ? &nlh
.nabnsb
: NULL
,
2101 terminate_now
= terminate
;
2105 wallcycle_stop(wcycle
,ewcMoveE
);
2106 bSumEkinhOld
= FALSE
;
2109 /* This is just for testing. Nothing is actually done to Ekin
2110 * since that would require extra communication.
2112 if (!bNEMD
&& debug
&& (vcm
->nr
> 0))
2115 mdatoms
->start
,mdatoms
->start
+mdatoms
->homenr
,
2116 state
->v
,vcm
->group_p
[0],
2117 mdatoms
->massT
,mdatoms
->tmass
,ekin
);
2120 /* Do center of mass motion removal */
2121 if (bStopCM
&& !bFFscan
&& !bRerunMD
)
2123 check_cm_grp(fplog
,vcm
,ir
,1);
2124 do_stopcm_grp(fplog
,mdatoms
->start
,mdatoms
->homenr
,mdatoms
->cVCM
,
2125 state
->x
,state
->v
,vcm
);
2126 inc_nrnb(nrnb
,eNR_STOPCM
,mdatoms
->homenr
);
2128 calc_vcm_grp(fplog,START(nsb),HOMENR(nsb),mdatoms->massT,x,v,vcm);
2129 check_cm_grp(fplog,vcm,ir);
2130 do_stopcm_grp(fplog,START(nsb),HOMENR(nsb),x,v,vcm);
2131 check_cm_grp(fplog,vcm,ir);
2135 /* Add force and shake contribution to the virial */
2136 m_add(force_vir
,shake_vir
,total_vir
);
2138 /* Calculate the amplitude of the cosine velocity profile */
2139 ekind
->cosacc
.vcos
= ekind
->cosacc
.mvcos
/mdatoms
->tmass
;
2141 /* Sum the kinetic energies of the groups & calc temp */
2142 enerd
->term
[F_TEMP
] = sum_ekin((bRerunMD
&& !rerun_fr
.bV
),
2143 &(ir
->opts
),ekind
,ekin
,
2144 &(enerd
->term
[F_DKDL
]));
2145 enerd
->term
[F_EKIN
] = trace(ekin
);
2147 /* Calculate pressure and apply LR correction if PPPM is used.
2148 * Use the box from last timestep since we already called update().
2150 enerd
->term
[F_PRES
] =
2151 calc_pres(fr
->ePBC
,ir
->nwall
,lastbox
,ekin
,total_vir
,pres
,
2152 (fr
->eeltype
==eelPPPM
)?enerd
->term
[F_COUL_RECIP
]:0.0);
2154 /* Calculate long range corrections to pressure and energy */
2155 if (bTCR
|| bFFscan
)
2157 set_avcsixtwelve(fplog
,fr
,top_global
);
2160 /* Calculate long range corrections to pressure and energy */
2161 calc_dispcorr(fplog
,ir
,fr
,step
,top_global
->natoms
,
2162 lastbox
,state
->lambda
,pres
,total_vir
,enerd
);
2164 sum_dhdl(enerd
,state
->lambda
,ir
);
2166 enerd
->term
[F_ETOT
] = enerd
->term
[F_EPOT
] + enerd
->term
[F_EKIN
];
2170 if (isnan(enerd
->term
[F_ETOT
]))
2171 gmx_fatal(FARGS
, "NaN detected at step %d\n",step
);
2174 if (_isnan(enerd
->term
[F_ETOT
]))
2175 gmx_fatal(FARGS
, "NaN detected at step %d\n",step
);
2184 enerd
->term
[F_ECONSERVED
] =
2185 enerd
->term
[F_ETOT
] +
2186 nosehoover_energy(&(ir
->opts
),ekind
,
2187 state
->nosehoover_xi
,
2188 state
->therm_integral
);
2191 enerd
->term
[F_ECONSERVED
] =
2192 enerd
->term
[F_ETOT
] +
2193 vrescale_energy(&(ir
->opts
),
2194 state
->therm_integral
);
2198 /* We can not just use bCalcEner, since then the simulation results
2199 * would depend on nstenergy and nstlog or step_nscheck.
2201 if ((state
->flags
& (1<<estPRES_PREV
)) && bNstEner
)
2203 /* Store the pressure in t_state for pressure coupling
2204 * at the next MD step.
2206 copy_mat(pres
,state
->pres_prev
);
2209 /* Check for excessively large energies */
2213 real etot_max
= 1e200
;
2215 real etot_max
= 1e30
;
2217 if (fabs(enerd
->term
[F_ETOT
]) > etot_max
)
2219 fprintf(stderr
,"Energy too large (%g), giving up\n",
2220 enerd
->term
[F_ETOT
]);
2226 /* The coordinates (x) were unshifted in update */
2227 if (bFFscan
&& (shellfc
==NULL
|| bConverged
))
2229 if (print_forcefield(fplog
,enerd
->term
,mdatoms
->homenr
,
2231 &(top_global
->mols
),mdatoms
->massT
,pres
))
2233 if (gmx_parallel_env())
2235 fprintf(stderr
,"\n");
2242 /* Only do GCT when the relaxation of shells (minimization) has converged,
2243 * otherwise we might be coupling to bogus energies.
2244 * In parallel we must always do this, because the other sims might
2248 /* Since this is called with the new coordinates state->x, I assume
2249 * we want the new box state->box too. / EL 20040121
2251 do_coupling(fplog
,oenv
,nfile
,fnm
,tcr
,t
,step
,enerd
->term
,fr
,
2253 mdatoms
,&(top
->idef
),mu_aver
,
2254 top_global
->mols
.nr
,cr
,
2255 state
->box
,total_vir
,pres
,
2256 mu_tot
,state
->x
,f
,bConverged
);
2260 /* Time for performance */
2261 if (((step
% stepout
) == 0) || bLastStep
)
2263 runtime_upd_proc(runtime
);
2273 upd_mdebin(mdebin
,bDoDHDL
? fp_dhdl
: NULL
,TRUE
,
2274 t
,mdatoms
->tmass
,enerd
,state
,lastbox
,
2275 shake_vir
,force_vir
,total_vir
,pres
,
2276 ekind
,mu_tot
,constr
);
2280 upd_mdebin_step(mdebin
);
2283 do_dr
= do_per_step(step
,ir
->nstdisreout
);
2284 do_or
= do_per_step(step
,ir
->nstorireout
);
2286 print_ebin(fp_ene
,do_ene
,do_dr
,do_or
,do_log
?fplog
:NULL
,step
,t
,
2287 eprNORMAL
,bCompact
,mdebin
,fcd
,groups
,&(ir
->opts
));
2289 if (ir
->ePull
!= epullNO
)
2291 pull_print_output(ir
->pull
,step
,t
);
2294 if (do_per_step(step
,ir
->nstlog
))
2296 if(fflush(fplog
) != 0)
2298 gmx_fatal(FARGS
,"Cannot flush logfile - maybe you are out of quota?");
2304 /* Remaining runtime */
2305 if (MULTIMASTER(cr
) && do_verbose
)
2309 fprintf(stderr
,"\n");
2311 print_time(stderr
,runtime
,step
,ir
);
2314 /* Replica exchange */
2316 if ((repl_ex_nst
> 0) && (step
> 0) && !bLastStep
&&
2317 do_per_step(step
,repl_ex_nst
))
2319 bExchanged
= replica_exchange(fplog
,cr
,repl_ex
,
2320 state_global
,enerd
->term
,
2323 if (bExchanged
&& PAR(cr
))
2325 if (DOMAINDECOMP(cr
))
2327 dd_partition_system(fplog
,step
,cr
,TRUE
,1,
2328 state_global
,top_global
,ir
,
2329 state
,&f
,mdatoms
,top
,fr
,
2330 vsite
,shellfc
,constr
,
2335 bcast_state(cr
,state
,FALSE
);
2343 /* read next frame from input trajectory */
2344 bNotLastFrame
= read_next_frame(oenv
,status
,&rerun_fr
);
2347 if (!bRerunMD
|| !rerun_fr
.bStep
)
2349 /* increase the MD step number */
2354 cycles
= wallcycle_stop(wcycle
,ewcSTEP
);
2355 if (DOMAINDECOMP(cr
) && wcycle
)
2357 dd_cycles_add(cr
->dd
,cycles
,ddCyclStep
);
2360 if (step_rel
== wcycle_get_reset_counters(wcycle
))
2362 /* Reset all the counters related to performance over the run */
2363 reset_all_counters(fplog
,cr
,step
,&step_rel
,ir
,wcycle
,nrnb
,runtime
);
2364 wcycle_set_reset_counters(wcycle
, 0);
2367 /* End of main MD loop */
2371 runtime_end(runtime
);
2378 if (!(cr
->duty
& DUTY_PME
))
2380 /* Tell the PME only node to finish */
2386 print_ebin(fp_ene
,FALSE
,FALSE
,FALSE
,fplog
,step
,t
,
2387 eprAVER
,FALSE
,mdebin
,fcd
,groups
,&(ir
->opts
));
2397 gmx_fio_fclose(fp_dhdl
);
2401 gmx_fio_fclose(fp_field
);
2406 if (ir
->nstlist
== -1 && nlh
.nns
> 0 && fplog
)
2408 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
)));
2409 fprintf(fplog
,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh
.ab
/nlh
.nns
);
2412 if (shellfc
&& fplog
)
2414 fprintf(fplog
,"Fraction of iterations that converged: %.2f %%\n",
2415 (nconverged
*100.0)/step_rel
);
2416 fprintf(fplog
,"Average number of force evaluations per MD step: %.2f\n\n",
2420 if (repl_ex_nst
> 0 && MASTER(cr
))
2422 print_replica_exchange_statistics(fplog
,repl_ex
);
2425 runtime
->nsteps_done
= step_rel
;