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"
95 // Temporary definitions - MRS //
96 #define CONVERGEITER 0.00000002
97 #define MAXITERCONST 500
100 #include "corewrap.h"
103 /* The following two variables and the signal_handler function
104 * is used from pme.c as well
106 extern bool bGotTermSignal
, bGotUsr1Signal
;
108 static RETSIGTYPE
signal_handler(int n
)
112 bGotTermSignal
= TRUE
;
116 bGotUsr1Signal
= TRUE
;
123 gmx_integrator_t
*func
;
126 /* The array should match the eI array in include/types/enums.h */
127 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
}};
129 /* Static variables for temporary use with the deform option */
130 static int init_step_tpx
;
131 static matrix box_tpx
;
133 static tMPI_Thread_mutex_t box_mutex
=TMPI_THREAD_MUTEX_INITIALIZER
;
136 bool done_iterating(bool *bFirstIterate
, bool *bIterate
, real fom
, real
*newf
, int n
)
139 /* monitor convergence, and use a secant search to propose new
142 The secant method computes x_{i+1} = x_{i} - f(x_{i}) * ---------------------
143 f(x_{i}) - f(x_{i-1})
145 The function we are trying to zero is fom-x, where fom is the
146 "figure of merit" which is the pressure (or the veta value) we
147 would get by putting in an old value of the pressure or veta into
148 the incrementor function for the step or half step. I have
149 verified that this gives the same answer as self consistent
150 iteration, usually in many fewer steps, especially for small tau_p.
152 We could possibly eliminate an iteration with proper use
153 of the value from the previous step, but that would take a bit
154 more bookkeeping, especially for veta, since tests indicate the
155 function of veta on the last step is not sufficiently close to
156 guarantee convergence this step. This is
157 good enough for now. On my tests, I could use tau_p down to
158 0.02, which is smaller that would ever be necessary in
159 practice. Generally, 3-5 iterations will be sufficient */
161 static real f
,fprev
,x
,xprev
;
167 *bFirstIterate
= FALSE
;
175 f
= fom
-x
; /* we want to zero this difference */
176 if ((iter_i
> 1) && (iter_i
< MAXITERCONST
))
184 *newf
= x
- (x
-xprev
)*(f
)/(f
-fprev
);
189 /* just use self-consistent iteration the first step to initialize, or
190 if it's not converging (which happens occasionally -- need to investigate why) */
194 /* Consider a slight shortcut allowing us to exit one sooner -- we check the
195 difference between the closest of x and xprev to the new
196 value. To be 100% certain, we should check the difference between
197 the last result, and the previous result, or
199 relerr = (fabs((x-xprev)/fom));
201 but this is pretty much never necessary under typical conditions.
202 Checking numerically, it seems to lead to almost exactly the same
203 trajectories, but there are small differences out a few decimal
204 places in the pressure, and eventually in the v_eta, but it could
207 if (fabs(*newf-x) < fabs(*newf - xprev)) { xmin = x;} else { xmin = xprev;}
208 relerr = (fabs((*newf-xmin) / *newf));
211 relerr
= (fabs((f
-fprev
)/fom
));
217 fprintf(debug
,"Iterating NPT constraints #%i: %6i %20.12f%14.6g%20.12f\n",n
,iter_i
,fom
,relerr
,*newf
);
220 if ((relerr
< CONVERGEITER
) || (fom
==0))
225 fprintf(debug
,"Iterating NPT constraints #%i: CONVERGED\n",n
);
229 if (iter_i
> 10*MAXITERCONST
)
231 gmx_fatal(FARGS
,"Could not converge NPT constraints\n");
243 void copy_coupling_state (t_state
*statea
,t_state
*stateb
,
244 tensor shake_vira
, tensor shake_virb
,
245 tensor ekina
, tensor ekinb
,
246 gmx_ekindata_t
*ekinda
,gmx_ekindata_t
*ekindb
)
249 /* MRS note -- might be able to get rid of some of the arguments. Look over it when it's all debugged */
253 if (ekina
) copy_mat(ekina
,ekinb
);
254 if (shake_vira
) copy_mat(shake_vira
,shake_virb
);
255 stateb
->natoms
= statea
->natoms
;
256 stateb
->ngtc
= statea
->ngtc
;
257 stateb
->veta
= statea
->veta
;
260 for (i
=0; i
<stateb
->ngtc
; i
++)
262 ekindb
->tcstat
[i
].T
= ekinda
->tcstat
[i
].T
;
263 ekindb
->tcstat
[i
].Th
= ekinda
->tcstat
[i
].Th
;
264 copy_mat(ekinda
->tcstat
[i
].ekinh
,ekindb
->tcstat
[i
].ekinh
);
265 copy_mat(ekinda
->tcstat
[i
].ekin
,ekindb
->tcstat
[i
].ekin
);
268 copy_rvecn(statea
->x
,stateb
->x
,0,stateb
->natoms
);
269 copy_rvecn(statea
->v
,stateb
->v
,0,stateb
->natoms
);
270 copy_mat(statea
->box
,stateb
->box
);
271 copy_mat(statea
->box_rel
,stateb
->box_rel
);
272 copy_mat(statea
->boxv
,stateb
->boxv
);
273 ngtc_eff
= stateb
->ngtc
+1;
274 /* need extra term for the barostat */
275 for (i
= 0; i
< ngtc_eff
; i
++)
277 for (j
=0; j
< NNHCHAIN
; j
++)
279 stateb
->nosehoover_xi
[i
*NNHCHAIN
+ j
] = statea
->nosehoover_xi
[i
*NNHCHAIN
+ j
];
280 stateb
->nosehoover_vxi
[i
*NNHCHAIN
+ j
] = statea
->nosehoover_vxi
[i
*NNHCHAIN
+ j
];
285 void compute_globals(FILE *fplog
, gmx_global_stat_t gstat
, t_commrec
*cr
, t_inputrec
*ir
,
286 t_forcerec
*fr
, gmx_ekindata_t
*ekind
,
287 t_state
*state
, t_state
*state_global
, t_mdatoms
*mdatoms
,
288 t_nrnb
*nrnb
, t_vcm
*vcm
, gmx_wallcycle_t wcycle
,
289 gmx_enerdata_t
*enerd
,tensor force_vir
, tensor shake_vir
,
290 tensor shakev_vir
, tensor total_vir
, tensor ekin
,
291 tensor pres
, rvec mu_tot
, gmx_constr_t constr
,
292 real
*chkpt
,real
*terminate
, real
*terminate_now
,
293 int *nabnsb
, matrix box
, gmx_mtop_t
*top_global
,real
*pcurr
,
294 bool *bSumEkinhOld
, bool bRerunMD
, bool bEkinFullStep
,
295 bool bDoCM
, bool bGStat
, bool bNEMD
, bool bFirstHalf
,
296 bool bIterate
, bool bFirstIterate
, bool bInitialize
,
297 bool bReadEkin
, int natoms
)
301 bool bTemp
=FALSE
,bPres
=FALSE
, bEner
=FALSE
;
302 real prescorr
,enercorr
,dvdlcorr
;
303 tensor corr_vir
,corr_pres
,shakeall_vir
;
305 /* decide when to calculate temperature and pressure.
306 Ener is only for calculating the dispersion correction. */
308 /* in initalization, it sums the shake virial in vv, and to
309 sums ekinh_old in leapfrog (or if we are calculating ekinh_old for other reasons */
330 if (bIterate
|| (ir
->epc
== epcTROTTER
))
356 /* ########## Kinetic energy ############## */
360 /* Non-equilibrium MD: this is parallellized, but only does communication
361 * when there really is NEMD.
364 if (PAR(cr
) && bNEMD
)
366 accumulate_u(cr
,&(ir
->opts
),ekind
);
371 restore_ekinstate_from_state(cr
,ekind
,&state_global
->ekinstate
);
375 calc_ke_part(state
,&(ir
->opts
),mdatoms
,ekind
,nrnb
,bEkinFullStep
);
379 /* Calculate center of mass velocity if necessary, also parallellized */
380 if (bDoCM
&& !bRerunMD
)
382 calc_vcm_grp(fplog
,mdatoms
->start
,mdatoms
->homenr
,mdatoms
,
383 state
->x
,state
->v
,vcm
);
391 GMX_MPE_LOG(ev_global_stat_start
);
393 global_stat(fplog
,gstat
,cr
,enerd
,force_vir
,shake_vir
,mu_tot
,
394 ir
,ekind
,FALSE
,bEkinFullStep
,constr
,vcm
,
395 NULL
,NULL
,terminate
,top_global
,state
);
396 GMX_MPE_LOG(ev_global_stat_finish
);
402 /* Sum the kinetic energies of the groups & calc temp */
403 enerd
->term
[F_TEMP
] = sum_ekin(bInitialize
,&(ir
->opts
),ekind
,ekin
,
404 &(enerd
->term
[F_DKDL
]),bEkinFullStep
);
405 enerd
->term
[F_EKIN
] = trace(ekin
);
408 /* ########## Long range energy data ###### */
412 calc_dispcorr(fplog
,ir
,fr
,0,top_global
,box
,state
->lambda
,
413 corr_pres
,corr_vir
,&prescorr
,&enercorr
,&dvdlcorr
);
416 if (bEner
&& bFirstIterate
)
418 enerd
->term
[F_DISPCORR
] = enercorr
;
419 enerd
->term
[F_EPOT
] += enercorr
;
420 enerd
->term
[F_DVDL
] += dvdlcorr
;
423 /* ########## Now pressure ############## */
428 /* Add force and shake contribution to the virial */
429 m_add(shake_vir
,shakev_vir
,shakeall_vir
);
430 m_add(force_vir
,shakeall_vir
,total_vir
);
432 /* Calculate pressure and apply LR correction if PPPM is used.
433 * Use the box from last timestep since we already called update().
436 enerd
->term
[F_PRES
] = calc_pres(fr
->ePBC
,ir
->nwall
,box
,ekin
,total_vir
,pres
,
437 (fr
->eeltype
==eelPPPM
)?enerd
->term
[F_COUL_RECIP
]:0.0);
439 /* Calculate long range corrections to pressure and energy */
440 /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
441 and computes enerd->term[F_DISPCORR]. Also modifies the
442 total_vir and pres tesors */
444 m_add(total_vir
,corr_vir
,total_vir
);
445 m_add(pres
,corr_pres
,pres
);
446 enerd
->term
[F_PDISPCORR
] = prescorr
;
447 enerd
->term
[F_PRES
] += prescorr
;
448 *pcurr
= enerd
->term
[F_PRES
];
454 struct mdrunner_arglist
468 const char *dddlb_opt
;
481 int ret
; /* return value */
485 static void mdrunner_start_fn(void *arg
)
487 struct mdrunner_arglist
*mda
=(struct mdrunner_arglist
*)arg
;
488 struct mdrunner_arglist mc
=*mda
; /* copy the arg list to make sure
489 that it's thread-local. This doesn't
490 copy pointed-to items, of course,
491 but those are all const. */
492 t_commrec
*cr
; /* we need a local version of this */
494 t_filenm
*fnm
=dup_tfn(mc
.nfile
, mc
.fnm
);
496 cr
=init_par_threads(mc
.cr
);
503 mda
->ret
=mdrunner(fplog
, cr
, mc
.nfile
, mc
.fnm
, mc
.oenv
, mc
.bVerbose
,
504 mc
.bCompact
, mc
.nstglobalcomm
,
505 mc
.ddxyz
, mc
.dd_node_order
, mc
.rdd
,
506 mc
.rconstr
, mc
.dddlb_opt
, mc
.dlb_scale
,
507 mc
.ddcsx
, mc
.ddcsy
, mc
.ddcsz
, mc
.nstepout
, mc
.nmultisim
,
508 mc
.repl_ex_nst
, mc
.repl_ex_seed
, mc
.pforce
,
509 mc
.cpt_period
, mc
.max_hours
, mc
.Flags
);
514 int mdrunner_threads(int nthreads
,
515 FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
516 const output_env_t oenv
, bool bVerbose
,bool bCompact
,
518 ivec ddxyz
,int dd_node_order
,real rdd
,real rconstr
,
519 const char *dddlb_opt
,real dlb_scale
,
520 const char *ddcsx
,const char *ddcsy
,const char *ddcsz
,
521 int nstepout
,int nmultisim
,int repl_ex_nst
,
522 int repl_ex_seed
, real pforce
,real cpt_period
,
523 real max_hours
, unsigned long Flags
)
526 /* first check whether we even need to start tMPI */
529 ret
=mdrunner(fplog
, cr
, nfile
, fnm
, oenv
, bVerbose
, bCompact
,
531 ddxyz
, dd_node_order
, rdd
, rconstr
, dddlb_opt
, dlb_scale
,
532 ddcsx
, ddcsy
, ddcsz
, nstepout
, nmultisim
, repl_ex_nst
,
533 repl_ex_seed
, pforce
, cpt_period
, max_hours
, Flags
);
538 struct mdrunner_arglist mda
;
539 /* fill the data structure to pass as void pointer to thread start fn */
545 mda
.bVerbose
=bVerbose
;
546 mda
.bCompact
=bCompact
;
547 mda
.nstglobalcomm
=nstglobalcomm
;
548 mda
.ddxyz
[XX
]=ddxyz
[XX
];
549 mda
.ddxyz
[YY
]=ddxyz
[YY
];
550 mda
.ddxyz
[ZZ
]=ddxyz
[ZZ
];
551 mda
.dd_node_order
=dd_node_order
;
554 mda
.dddlb_opt
=dddlb_opt
;
555 mda
.dlb_scale
=dlb_scale
;
559 mda
.nstepout
=nstepout
;
560 mda
.nmultisim
=nmultisim
;
561 mda
.repl_ex_nst
=repl_ex_nst
;
562 mda
.repl_ex_seed
=repl_ex_seed
;
564 mda
.cpt_period
=cpt_period
;
565 mda
.max_hours
=max_hours
;
568 fprintf(stderr
, "Starting %d threads\n",nthreads
);
570 tMPI_Init_fn(nthreads
, mdrunner_start_fn
, (void*)(&mda
) );
574 gmx_comm("Multiple threads requested but not compiled with threads");
581 int mdrunner(FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
582 const output_env_t oenv
, bool bVerbose
,bool bCompact
,
584 ivec ddxyz
,int dd_node_order
,real rdd
,real rconstr
,
585 const char *dddlb_opt
,real dlb_scale
,
586 const char *ddcsx
,const char *ddcsy
,const char *ddcsz
,
587 int nstepout
,int nmultisim
,int repl_ex_nst
,int repl_ex_seed
,
588 real pforce
,real cpt_period
,real max_hours
,
591 double nodetime
=0,realtime
;
592 t_inputrec
*inputrec
;
599 gmx_mtop_t
*mtop
=NULL
;
600 t_mdatoms
*mdatoms
=NULL
;
604 gmx_pme_t
*pmedata
=NULL
;
605 gmx_vsite_t
*vsite
=NULL
;
607 int i
,m
,nChargePerturbed
=-1,status
,nalloc
;
609 gmx_wallcycle_t wcycle
;
610 bool bReadRNG
,bReadEkin
;
612 gmx_runtime_t runtime
;
614 gmx_large_int_t reset_counters
;
617 /* Essential dynamics */
618 if (opt2bSet("-ei",nfile
,fnm
))
620 /* Open input and output files, allocate space for ED data structure */
621 ed
= ed_open(nfile
,fnm
,cr
);
629 if (bVerbose
&& SIMMASTER(cr
))
631 fprintf(stderr
,"Getting Loaded...\n");
634 if (Flags
& MD_APPENDFILES
)
641 /* The master thread on the master node reads from disk,
642 * then distributes everything to the other processors.
645 list
= (SIMMASTER(cr
) && !(Flags
& MD_APPENDFILES
)) ? (LIST_SCALARS
| LIST_INPUTREC
) : 0;
648 init_parallel(fplog
, opt2fn_master("-s",nfile
,fnm
,cr
),cr
,
649 inputrec
,mtop
,state
,list
);
654 /* Read a file for a single processor */
656 init_single(fplog
,inputrec
,ftp2fn(efTPX
,nfile
,fnm
),mtop
,state
);
658 if (!EEL_PME(inputrec
->coulombtype
) || (Flags
& MD_PARTDEC
))
664 fcRegisterSteps(inputrec
->nsteps
,inputrec
->init_step
);
667 /* NMR restraints must be initialized before load_checkpoint,
668 * since with time averaging the history is added to t_state.
669 * For proper consistency check we therefore need to extend
671 * So the PME-only nodes (if present) will also initialize
672 * the distance restraints.
676 /* This needs to be called before read_checkpoint to extend the state */
677 init_disres(fplog
,mtop
,inputrec
,cr
,Flags
& MD_PARTDEC
,fcd
,state
);
679 if (gmx_mtop_ftype_count(mtop
,F_ORIRES
) > 0)
681 if (PAR(cr
) && !(Flags
& MD_PARTDEC
))
683 gmx_fatal(FARGS
,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
685 /* Orientation restraints */
688 init_orires(fplog
,mtop
,state
->x
,inputrec
,cr
->ms
,&(fcd
->orires
),
693 if (DEFORM(*inputrec
))
695 /* Store the deform reference box before reading the checkpoint */
698 copy_mat(state
->box
,box
);
702 gmx_bcast(sizeof(box
),box
,cr
);
704 /* Because we do not have the update struct available yet
705 * in which the reference values should be stored,
706 * we store them temporarily in static variables.
707 * This should be thread safe, since they are only written once
708 * and with identical values.
711 tMPI_Thread_mutex_lock(&box_mutex
);
713 init_step_tpx
= inputrec
->init_step
;
714 copy_mat(box
,box_tpx
);
716 tMPI_Thread_mutex_unlock(&box_mutex
);
720 if (opt2bSet("-cpi",nfile
,fnm
))
722 /* Check if checkpoint file exists before doing continuation.
723 * This way we can use identical input options for the first and subsequent runs...
725 if( gmx_fexist_master(opt2fn_master("-cpi",nfile
,fnm
,cr
),cr
) )
727 load_checkpoint(opt2fn_master("-cpi",nfile
,fnm
,cr
),fplog
,
728 cr
,Flags
& MD_PARTDEC
,ddxyz
,
729 inputrec
,state
,&bReadRNG
,&bReadEkin
,
730 (Flags
& MD_APPENDFILES
));
734 Flags
|= MD_READ_RNG
;
738 Flags
|= MD_READ_EKIN
;
743 if (MASTER(cr
) && (Flags
& MD_APPENDFILES
))
745 fplog
= gmx_log_open(ftp2fn(efLOG
,nfile
,fnm
),cr
,!(Flags
& MD_SEPPOT
),
751 copy_mat(state
->box
,box
);
756 gmx_bcast(sizeof(box
),box
,cr
);
759 if (bVerbose
&& SIMMASTER(cr
))
761 fprintf(stderr
,"Loaded with Money\n\n");
764 if (PAR(cr
) && !((Flags
& MD_PARTDEC
) || EI_TPI(inputrec
->eI
)))
766 cr
->dd
= init_domain_decomposition(fplog
,cr
,Flags
,ddxyz
,rdd
,rconstr
,
773 make_dd_communicators(fplog
,cr
,dd_node_order
);
775 /* Set overallocation to avoid frequent reallocation of arrays */
776 set_over_alloc_dd(TRUE
);
780 cr
->duty
= (DUTY_PP
| DUTY_PME
);
781 npme_major
= cr
->nnodes
;
783 if (inputrec
->ePBC
== epbcSCREW
)
786 "pbc=%s is only implemented with domain decomposition",
787 epbc_names
[inputrec
->ePBC
]);
793 /* After possible communicator splitting in make_dd_communicators.
794 * we can set up the intra/inter node communication.
796 gmx_setup_nodecomm(fplog
,cr
);
799 wcycle
= wallcycle_init(fplog
,cr
);
802 /* Master synchronizes its value of reset_counters with all nodes
803 * including PME only nodes */
804 reset_counters
= wcycle_get_reset_counters(wcycle
);
805 gmx_bcast_sim(sizeof(reset_counters
),&reset_counters
,cr
);
806 wcycle_set_reset_counters(wcycle
, reset_counters
);
811 if (cr
->duty
& DUTY_PP
)
813 /* For domain decomposition we allocate dynamically
814 * in dd_partition_system.
816 if (DOMAINDECOMP(cr
))
818 bcast_state_setup(cr
,state
);
828 bcast_state(cr
,state
,TRUE
);
832 /* Dihedral Restraints */
833 if (gmx_mtop_ftype_count(mtop
,F_DIHRES
) > 0)
835 init_dihres(fplog
,mtop
,inputrec
,fcd
);
838 /* Initiate forcerecord */
840 init_forcerec(fplog
,oenv
,fr
,fcd
,inputrec
,mtop
,cr
,box
,FALSE
,
841 opt2fn("-table",nfile
,fnm
),
842 opt2fn("-tablep",nfile
,fnm
),
843 opt2fn("-tableb",nfile
,fnm
),FALSE
,pforce
);
845 /* version for PCA_NOT_READ_NODE (see md.c) */
846 /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
847 "nofile","nofile","nofile",FALSE,pforce);
849 fr
->bSepDVDL
= ((Flags
& MD_SEPPOT
) == MD_SEPPOT
);
851 /* Initialize QM-MM */
854 init_QMMMrec(cr
,box
,mtop
,inputrec
,fr
);
857 /* Initialize the mdatoms structure.
858 * mdatoms is not filled with atom data,
859 * as this can not be done now with domain decomposition.
861 mdatoms
= init_mdatoms(fplog
,mtop
,inputrec
->efep
!=efepNO
);
863 /* Initialize the virtual site communication */
864 vsite
= init_vsite(mtop
,cr
);
866 calc_shifts(box
,fr
->shift_vec
);
868 /* With periodic molecules the charge groups should be whole at start up
869 * and the virtual sites should not be far from their proper positions.
871 if (!inputrec
->bContinuation
&& MASTER(cr
) &&
872 !(inputrec
->ePBC
!= epbcNONE
&& inputrec
->bPeriodicMols
))
874 /* Make molecules whole at start of run */
875 if (fr
->ePBC
!= epbcNONE
)
877 do_pbc_first_mtop(fplog
,inputrec
->ePBC
,box
,mtop
,state
->x
);
881 /* Correct initial vsite positions are required
882 * for the initial distribution in the domain decomposition
883 * and for the initial shell prediction.
885 construct_vsites_mtop(fplog
,vsite
,mtop
,state
->x
);
889 /* Initiate PPPM if necessary */
890 if (fr
->eeltype
== eelPPPM
)
892 if (mdatoms
->nChargePerturbed
)
894 gmx_fatal(FARGS
,"Free energy with %s is not implemented",
895 eel_names
[fr
->eeltype
]);
897 status
= gmx_pppm_init(fplog
,cr
,oenv
,FALSE
,TRUE
,box
,
898 getenv("GMXGHAT"),inputrec
, (Flags
& MD_REPRODUCIBLE
));
901 gmx_fatal(FARGS
,"Error %d initializing PPPM",status
);
905 if (EEL_PME(fr
->eeltype
))
907 ewaldcoeff
= fr
->ewaldcoeff
;
908 pmedata
= &fr
->pmedata
;
917 /* This is a PME only node */
919 /* We don't need the state */
922 ewaldcoeff
= calc_ewaldcoeff(inputrec
->rcoulomb
, inputrec
->ewald_rtol
);
926 /* Initiate PME if necessary,
927 * either on all nodes or on dedicated PME nodes only. */
928 if (EEL_PME(inputrec
->coulombtype
))
932 nChargePerturbed
= mdatoms
->nChargePerturbed
;
934 if (cr
->npmenodes
> 0)
936 /* The PME only nodes need to know nChargePerturbed */
937 gmx_bcast_sim(sizeof(nChargePerturbed
),&nChargePerturbed
,cr
);
939 if (cr
->duty
& DUTY_PME
)
941 status
= gmx_pme_init(pmedata
,cr
,npme_major
,inputrec
,
942 mtop
? mtop
->natoms
: 0,nChargePerturbed
,
943 (Flags
& MD_REPRODUCIBLE
));
946 gmx_fatal(FARGS
,"Error %d initializing PME",status
);
952 if (integrator
[inputrec
->eI
].func
== do_md
)
954 /* Turn on signal handling on all nodes */
956 * (A user signal from the PME nodes (if any)
957 * is communicated to the PP nodes.
959 if (getenv("GMX_NO_TERM") == NULL
)
963 fprintf(debug
,"Installing signal handler for SIGTERM\n");
965 signal(SIGTERM
,signal_handler
);
968 if (getenv("GMX_NO_USR1") == NULL
)
972 fprintf(debug
,"Installing signal handler for SIGUSR1\n");
974 signal(SIGUSR1
,signal_handler
);
979 if (cr
->duty
& DUTY_PP
)
981 if (inputrec
->ePull
!= epullNO
)
983 /* Initialize pull code */
984 init_pull(fplog
,inputrec
,nfile
,fnm
,mtop
,cr
,oenv
,
985 EI_DYNAMICS(inputrec
->eI
) && MASTER(cr
),Flags
);
988 constr
= init_constraints(fplog
,mtop
,inputrec
,ed
,state
,cr
);
990 if (DOMAINDECOMP(cr
))
992 dd_init_bondeds(fplog
,cr
->dd
,mtop
,vsite
,constr
,inputrec
,
993 Flags
& MD_DDBONDCHECK
,fr
->cginfo_mb
);
995 set_dd_parameters(fplog
,cr
->dd
,dlb_scale
,inputrec
,fr
,&ddbox
);
997 setup_dd_grid(fplog
,cr
->dd
);
1000 /* Now do whatever the user wants us to do (how flexible...) */
1001 integrator
[inputrec
->eI
].func(fplog
,cr
,nfile
,fnm
,
1002 oenv
,bVerbose
,bCompact
,
1005 nstepout
,inputrec
,mtop
,
1007 mdatoms
,nrnb
,wcycle
,ed
,fr
,
1008 repl_ex_nst
,repl_ex_seed
,
1009 cpt_period
,max_hours
,
1013 if (inputrec
->ePull
!= epullNO
)
1015 finish_pull(fplog
,inputrec
->pull
);
1021 gmx_pmeonly(*pmedata
,cr
,nrnb
,wcycle
,ewaldcoeff
,FALSE
,inputrec
);
1024 if (EI_DYNAMICS(inputrec
->eI
) || EI_TPI(inputrec
->eI
))
1026 /* Some timing stats */
1029 if (runtime
.proc
== 0)
1031 runtime
.proc
= runtime
.real
;
1040 wallcycle_stop(wcycle
,ewcRUN
);
1042 /* Finish up, write some stuff
1043 * if rerunMD, don't write last frame again
1045 finish_run(fplog
,cr
,ftp2fn(efSTO
,nfile
,fnm
),
1046 inputrec
,nrnb
,wcycle
,&runtime
,
1047 EI_DYNAMICS(inputrec
->eI
) && !MULTISIM(cr
));
1049 /* Does what it says */
1050 print_date_and_time(fplog
,cr
->nodeid
,"Finished mdrun",&runtime
);
1052 /* Close logfile already here if we were appending to it */
1053 if (MASTER(cr
) && (Flags
& MD_APPENDFILES
))
1055 gmx_log_close(fplog
);
1062 else if(bGotUsr1Signal
)
1074 static void md_print_warning(const t_commrec
*cr
,FILE *fplog
,const char *buf
)
1078 fprintf(stderr
,"\n%s\n",buf
);
1082 fprintf(fplog
,"\n%s\n",buf
);
1086 static void check_nst_param(FILE *fplog
,t_commrec
*cr
,
1087 const char *desc_nst
,int nst
,
1088 const char *desc_p
,int *p
)
1092 if (*p
> 0 && *p
% nst
!= 0)
1094 /* Round up to the next multiple of nst */
1095 *p
= ((*p
)/nst
+ 1)*nst
;
1096 sprintf(buf
,"NOTE: %s changes %s to %d\n",desc_nst
,desc_p
,*p
);
1097 md_print_warning(cr
,fplog
,buf
);
1101 static void reset_all_counters(FILE *fplog
,t_commrec
*cr
,
1102 gmx_large_int_t step
,
1103 gmx_large_int_t
*step_rel
,t_inputrec
*ir
,
1104 gmx_wallcycle_t wcycle
,t_nrnb
*nrnb
,
1105 gmx_runtime_t
*runtime
)
1107 char buf
[STRLEN
],sbuf
[22];
1109 /* Reset all the counters related to performance over the run */
1110 sprintf(buf
,"Step %s: resetting all time and cycle counters\n",
1111 gmx_step_str(step
,sbuf
));
1112 md_print_warning(cr
,fplog
,buf
);
1114 wallcycle_stop(wcycle
,ewcRUN
);
1115 wallcycle_reset_all(wcycle
);
1116 if (DOMAINDECOMP(cr
))
1118 reset_dd_statistics_counters(cr
->dd
);
1121 ir
->init_step
+= *step_rel
;
1122 ir
->nsteps
-= *step_rel
;
1124 wallcycle_start(wcycle
,ewcRUN
);
1125 runtime_start(runtime
);
1126 print_date_and_time(fplog
,cr
->nodeid
,"Restarted time",runtime
);
1129 static int check_nstglobalcomm(FILE *fplog
,t_commrec
*cr
,
1130 int nstglobalcomm
,t_inputrec
*ir
)
1134 if (!EI_DYNAMICS(ir
->eI
))
1139 if (nstglobalcomm
== -1)
1141 if (ir
->nstcalcenergy
== 0 && ir
->nstlist
== 0)
1144 if (ir
->nstenergy
> 0 && ir
->nstenergy
< nstglobalcomm
)
1146 nstglobalcomm
= ir
->nstenergy
;
1151 /* We assume that if nstcalcenergy > nstlist,
1152 * nstcalcenergy is a multiple of nstlist.
1154 if (ir
->nstcalcenergy
== 0 ||
1155 (ir
->nstlist
> 0 && ir
->nstlist
< ir
->nstcalcenergy
))
1157 nstglobalcomm
= ir
->nstlist
;
1161 nstglobalcomm
= ir
->nstcalcenergy
;
1167 if (ir
->nstlist
> 0 &&
1168 nstglobalcomm
> ir
->nstlist
&& nstglobalcomm
% ir
->nstlist
!= 0)
1170 nstglobalcomm
= (nstglobalcomm
/ ir
->nstlist
)*ir
->nstlist
;
1171 sprintf(buf
,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm
);
1172 md_print_warning(cr
,fplog
,buf
);
1174 if (nstglobalcomm
> ir
->nstcalcenergy
)
1176 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
1177 "nstcalcenergy",&ir
->nstcalcenergy
);
1180 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
1181 "nstenergy",&ir
->nstenergy
);
1183 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
1184 "nstlog",&ir
->nstlog
);
1187 if (ir
->comm_mode
!= ecmNO
&& ir
->nstcomm
< nstglobalcomm
)
1189 sprintf(buf
,"WARNING: Changing nstcomm from %d to %d\n",
1190 ir
->nstcomm
,nstglobalcomm
);
1191 md_print_warning(cr
,fplog
,buf
);
1192 ir
->nstcomm
= nstglobalcomm
;
1195 return nstglobalcomm
;
1198 static void check_ir_old_tpx_versions(t_commrec
*cr
,FILE *fplog
,
1199 t_inputrec
*ir
,gmx_mtop_t
*mtop
)
1201 /* Check required for old tpx files */
1202 if (IR_TWINRANGE(*ir
) && ir
->nstlist
> 1 &&
1203 ir
->nstcalcenergy
% ir
->nstlist
!= 0)
1205 md_print_warning(cr
,fplog
,"Old tpr file with twin-range settings: modifiying energy calculation and/or T/P-coupling frequencies");
1207 if (gmx_mtop_ftype_count(mtop
,F_CONSTR
) +
1208 gmx_mtop_ftype_count(mtop
,F_CONSTRNC
) > 0 &&
1209 ir
->eConstrAlg
== econtSHAKE
)
1211 md_print_warning(cr
,fplog
,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect");
1212 if (ir
->epc
!= epcNO
)
1214 gmx_fatal(FARGS
,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
1217 check_nst_param(fplog
,cr
,"nstlist",ir
->nstlist
,
1218 "nstcalcenergy",&ir
->nstcalcenergy
);
1220 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
1221 "nstenergy",&ir
->nstenergy
);
1222 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
1223 "nstlog",&ir
->nstlog
);
1224 if (ir
->efep
!= efepNO
)
1226 check_nst_param(fplog
,cr
,"nstdhdl",ir
->nstdhdl
,
1227 "nstenergy",&ir
->nstenergy
);
1232 bool bGStatEveryStep
;
1233 gmx_large_int_t step_ns
;
1234 gmx_large_int_t step_nscheck
;
1235 gmx_large_int_t nns
;
1245 static void reset_nlistheuristics(gmx_nlheur_t
*nlh
,gmx_large_int_t step
)
1249 nlh
->step_nscheck
= step
;
1252 static void init_nlistheuristics(gmx_nlheur_t
*nlh
,
1253 bool bGStatEveryStep
,gmx_large_int_t step
)
1255 nlh
->bGStatEveryStep
= bGStatEveryStep
;
1262 reset_nlistheuristics(nlh
,step
);
1265 static void update_nliststatistics(gmx_nlheur_t
*nlh
,gmx_large_int_t step
)
1267 gmx_large_int_t nl_lt
;
1268 char sbuf
[22],sbuf2
[22];
1270 /* Determine the neighbor list life time */
1271 nl_lt
= step
- nlh
->step_ns
;
1274 fprintf(debug
,"%d atoms beyond ns buffer, updating neighbor list after %s steps\n",nlh
->nabnsb
,gmx_step_str(nl_lt
,sbuf
));
1278 nlh
->s2
+= nl_lt
*nl_lt
;
1279 nlh
->ab
+= nlh
->nabnsb
;
1280 if (nlh
->lt_runav
== 0)
1282 nlh
->lt_runav
= nl_lt
;
1283 /* Initialize the fluctuation average
1284 * such that at startup we check after 0 steps.
1286 nlh
->lt_runav2
= sqr(nl_lt
/2.0);
1288 /* Running average with 0.9 gives an exp. history of 9.5 */
1289 nlh
->lt_runav2
= 0.9*nlh
->lt_runav2
+ 0.1*sqr(nlh
->lt_runav
- nl_lt
);
1290 nlh
->lt_runav
= 0.9*nlh
->lt_runav
+ 0.1*nl_lt
;
1291 if (nlh
->bGStatEveryStep
)
1293 /* Always check the nlist validity */
1294 nlh
->step_nscheck
= step
;
1298 /* We check after: <life time> - 2*sigma
1299 * The factor 2 is quite conservative,
1300 * but we assume that with nstlist=-1 the user
1301 * prefers exact integration over performance.
1303 nlh
->step_nscheck
= step
1304 + (int)(nlh
->lt_runav
- 2.0*sqrt(nlh
->lt_runav2
)) - 1;
1308 fprintf(debug
,"nlist life time %s run av. %4.1f sig %3.1f check %s check with -gcom %d\n",
1309 gmx_step_str(nl_lt
,sbuf
),nlh
->lt_runav
,sqrt(nlh
->lt_runav2
),
1310 gmx_step_str(nlh
->step_nscheck
-step
+1,sbuf2
),
1311 (int)(nlh
->lt_runav
- 2.0*sqrt(nlh
->lt_runav2
)));
1315 static void set_nlistheuristics(gmx_nlheur_t
*nlh
,bool bReset
,gmx_large_int_t step
)
1321 reset_nlistheuristics(nlh
,step
);
1325 update_nliststatistics(nlh
,step
);
1328 nlh
->step_ns
= step
;
1329 /* Initialize the cumulative coordinate scaling matrix */
1330 clear_mat(nlh
->scale_tot
);
1331 for(d
=0; d
<DIM
; d
++)
1333 nlh
->scale_tot
[d
][d
] = 1.0;
1337 double do_md(FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
1338 const output_env_t oenv
, bool bVerbose
,bool bCompact
,
1340 gmx_vsite_t
*vsite
,gmx_constr_t constr
,
1341 int stepout
,t_inputrec
*ir
,
1342 gmx_mtop_t
*top_global
,
1344 t_state
*state_global
,
1346 t_nrnb
*nrnb
,gmx_wallcycle_t wcycle
,
1347 gmx_edsam_t ed
,t_forcerec
*fr
,
1348 int repl_ex_nst
,int repl_ex_seed
,
1349 real cpt_period
,real max_hours
,
1350 unsigned long Flags
,
1351 gmx_runtime_t
*runtime
)
1353 int fp_trn
=0,fp_xtc
=0;
1354 ener_file_t fp_ene
=NULL
;
1355 gmx_large_int_t step
,step_rel
;
1357 FILE *fp_dhdl
=NULL
,*fp_field
=NULL
;
1360 bool bGStatEveryStep
,bGStat
,bNstEner
,bCalcPres
,bCalcEner
;
1361 bool bNS
,bNStList
,bSimAnn
,bStopCM
,bRerunMD
,bNotLastFrame
=FALSE
,
1362 bFirstStep
,bStateFromTPX
,bInitStep
,bLastStep
,bEkinFullStep
,bBornRadii
;
1364 bool bNEMD
,do_ene
,do_log
,do_verbose
,bRerunWarnNoV
=TRUE
,
1365 bForceUpdate
=FALSE
,bX
,bV
,bF
,bXTC
,bCPT
;
1368 tensor force_vir
,shake_vir
,shake_Vir_save
,total_vir
,total_vir_new
,prev_vir
,pres
,ekin
,ekin_save
;
1373 matrix
*scale_tot
,pcoupl_mu
,M
,ebox
;
1375 t_trxframe rerun_fr
;
1376 gmx_repl_ex_t repl_ex
=NULL
;
1378 /* Booleans (disguised as a reals) to checkpoint and terminate mdrun */
1379 real chkpt
=0,terminate
=0,terminate_now
=0;
1381 gmx_localtop_t
*top
;
1382 t_mdebin
*mdebin
=NULL
;
1383 t_state
*state
=NULL
;
1384 rvec
*f_global
=NULL
;
1387 gmx_enerdata_t
*enerd
;
1389 gmx_global_stat_t gstat
;
1390 gmx_update_t upd
=NULL
;
1391 t_graph
*graph
=NULL
;
1394 gmx_groups_t
*groups
;
1395 gmx_ekindata_t
*ekind
, *ekind_save
;
1396 gmx_shellfc_t shellfc
;
1397 int count
,nconverged
=0;
1401 bool bTCR
=FALSE
,bConverged
=TRUE
,bOK
,bSumEkinhOld
,bExchanged
;
1403 book bIterate
,bFirstIterate
,bTemp
,bPres
,bTrotter
;
1404 real temp0
,mu_aver
=0,dvdl
;
1406 atom_id
*grpindex
=NULL
;
1408 t_coupl_rec
*tcr
=NULL
;
1409 rvec
*xcopy
=NULL
,*vcopy
=NULL
,*xbuf
=NULL
;
1410 matrix boxcopy
,lastbox
;
1412 real fom
,oldfom
,veta_save
,pcurr
,scalevir
,tracevir
;
1415 int reset_counters
=-1;
1416 real last_conserved
= 0;
1420 char sbuf
[22],sbuf2
[22];
1421 bool bHandledSignal
=FALSE
;
1423 /* Temporary addition for FAHCORE checkpointing */
1428 /* Check for special mdrun options */
1429 bRerunMD
= (Flags
& MD_RERUN
);
1430 bIonize
= (Flags
& MD_IONIZE
);
1431 bFFscan
= (Flags
& MD_FFSCAN
);
1432 bAppend
= (Flags
& MD_APPENDFILES
);
1433 /* all the iteratative cases - really, only if there are constraints */
1434 bIterate
= ((ir
->epc
== epcTROTTER
) && (constr
));
1435 bEkinFullStep
= ((ir
->eI
== eiVV
) && !(ir
->etc
==etcTROTTEREKINH
));
1436 bTrotter
= ((ir
->eI
== eiVV
) && ((ir
->etc
== etcTROTTER
) ||
1437 (ir
->etc
== etcTROTTEREKINH
) || (ir
->etc
== epcTROTTER
)));
1442 /* Since we don't know if the frames read are related in any way,
1443 * rebuild the neighborlist at every step.
1446 ir
->nstcalcenergy
= 1;
1450 check_ir_old_tpx_versions(cr
,fplog
,ir
,top_global
);
1452 nstglobalcomm
= check_nstglobalcomm(fplog
,cr
,nstglobalcomm
,ir
);
1453 bGStatEveryStep
= (nstglobalcomm
== 1);
1455 if (!bGStatEveryStep
&& ir
->nstlist
== -1 && fplog
!= NULL
)
1458 "To reduce the energy communication with nstlist = -1\n"
1459 "the neighbor list validity should not be checked at every step,\n"
1460 "this means that exact integration is not guaranteed.\n"
1461 "The neighbor list validity is checked after:\n"
1462 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
1463 "In most cases this will result in exact integration.\n"
1464 "This reduces the energy communication by a factor of 2 to 3.\n"
1465 "If you want less energy communication, set nstlist > 3.\n\n");
1468 if (bRerunMD
|| bFFscan
)
1472 groups
= &top_global
->groups
;
1474 /* Initial values */
1475 init_md(fplog
,cr
,ir
,oenv
,&t
,&t0
,&state_global
->lambda
,&lam0
,
1476 nrnb
,top_global
,&upd
,
1477 nfile
,fnm
,&fp_trn
,&fp_xtc
,&fp_ene
,&fn_cpt
,
1478 &fp_dhdl
,&fp_field
,&mdebin
,
1479 force_vir
,shake_vir
,mu_tot
,&bNEMD
,&bSimAnn
,&vcm
,Flags
);
1481 /* Energy terms and groups */
1483 clear_mat(shakev_vir
);
1484 init_enerdata(top_global
->groups
.grps
[egcENER
].nr
,ir
->n_flambda
,enerd
);
1485 if (DOMAINDECOMP(cr
))
1491 snew(f
,top_global
->natoms
);
1494 /* Kinetic energy data */
1496 init_ekindata(fplog
,top_global
,&(ir
->opts
),ekind
);
1497 /* needed for iteration of constraints */
1499 init_ekindata(fplog
,top_global
,&(ir
->opts
),ekind_save
);
1500 /* Copy the cos acceleration to the groups struct */
1501 ekind
->cosacc
.cos_accel
= ir
->cos_accel
;
1503 gstat
= global_stat_init(ir
);
1506 /* Check for polarizable models and flexible constraints */
1507 shellfc
= init_shell_flexcon(fplog
,
1508 top_global
,n_flexible_constraints(constr
),
1509 (ir
->bContinuation
||
1510 (DOMAINDECOMP(cr
) && !MASTER(cr
))) ?
1511 NULL
: state_global
->x
);
1516 tMPI_Thread_mutex_lock(&box_mutex
);
1518 set_deform_reference_box(upd
,init_step_tpx
,box_tpx
);
1520 tMPI_Thread_mutex_unlock(&box_mutex
);
1525 double io
= compute_io(ir
,top_global
->natoms
,groups
,mdebin
->ebin
->nener
,1);
1526 if ((io
> 2000) && MASTER(cr
))
1528 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
1532 if (DOMAINDECOMP(cr
)) {
1533 top
= dd_init_local_top(top_global
);
1536 dd_init_local_state(cr
->dd
,state_global
,state
);
1538 if (DDMASTER(cr
->dd
) && ir
->nstfout
) {
1539 snew(f_global
,state_global
->natoms
);
1543 /* Initialize the particle decomposition and split the topology */
1544 top
= split_system(fplog
,top_global
,ir
,cr
);
1546 pd_cg_range(cr
,&fr
->cg0
,&fr
->hcg
);
1547 pd_at_range(cr
,&a0
,&a1
);
1549 top
= gmx_mtop_generate_local_top(top_global
,ir
);
1552 a1
= top_global
->natoms
;
1555 state
= partdec_init_local_state(cr
,state_global
);
1558 atoms2md(top_global
,ir
,0,NULL
,a0
,a1
-a0
,mdatoms
);
1561 set_vsite_top(vsite
,top
,mdatoms
,cr
);
1564 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
) {
1565 graph
= mk_graph(fplog
,&(top
->idef
),0,top_global
->natoms
,FALSE
,FALSE
);
1569 make_local_shells(cr
,mdatoms
,shellfc
);
1572 if (ir
->pull
&& PAR(cr
)) {
1573 dd_make_local_pull_groups(NULL
,ir
->pull
,mdatoms
);
1577 if (DOMAINDECOMP(cr
))
1579 /* Distribute the charge groups over the nodes from the master node */
1580 dd_partition_system(fplog
,ir
->init_step
,cr
,TRUE
,1,
1581 state_global
,top_global
,ir
,
1582 state
,&f
,mdatoms
,top
,fr
,
1583 vsite
,shellfc
,constr
,
1587 /* If not DD, copy gb data */
1588 if(ir
->implicit_solvent
&& !DOMAINDECOMP(cr
))
1590 make_local_gb(cr
,fr
->born
,ir
->gb_algorithm
);
1593 update_mdatoms(mdatoms
,state
->lambda
);
1597 /* Update mdebin with energy history if appending to output files */
1598 if ( Flags
& MD_APPENDFILES
)
1600 restore_energyhistory_from_state(mdebin
,&state_global
->enerhist
);
1602 /* Set the initial energy history in state to zero by updating once */
1603 update_energyhistory(&state_global
->enerhist
,mdebin
);
1607 if ((state
->flags
& (1<<estLD_RNG
)) && (Flags
& MD_READ_RNG
)) {
1608 /* Set the random state if we read a checkpoint file */
1609 set_stochd_state(upd
,state
);
1612 /* Initialize constraints */
1614 if (!DOMAINDECOMP(cr
))
1615 set_constraints(constr
,top
,ir
,mdatoms
,cr
);
1618 /* Check whether we have to GCT stuff */
1619 bTCR
= ftp2bSet(efGCT
,nfile
,fnm
);
1622 fprintf(stderr
,"Will do General Coupling Theory!\n");
1624 gnx
= top_global
->mols
.nr
;
1626 for(i
=0; (i
<gnx
); i
++) {
1631 if (repl_ex_nst
> 0 && MASTER(cr
))
1632 repl_ex
= init_replica_exchange(fplog
,cr
->ms
,state_global
,ir
,
1633 repl_ex_nst
,repl_ex_seed
);
1636 if (!ir
->bContinuation
&& !bRerunMD
)
1638 if (mdatoms
->cFREEZE
&& (state
->flags
& (1<<estV
)))
1640 /* Set the velocities of frozen particles to zero */
1641 for(i
=mdatoms
->start
; i
<mdatoms
->start
+mdatoms
->homenr
; i
++)
1643 for(m
=0; m
<DIM
; m
++)
1645 if (ir
->opts
.nFreeze
[mdatoms
->cFREEZE
[i
]][m
])
1655 /* Constrain the initial coordinates and velocities */
1656 /* added calls to force to back-construct the constraint virial */
1657 force_flags
= (GMX_FORCE_STATECHANGED
|
1658 GMX_FORCE_ALLFORCES
|
1663 do_constrain_first(fplog
,constr
,ir
,mdatoms
,state
,f
,
1664 graph
,cr
,nrnb
,fr
,top
,top_global
,
1665 fcd
,wcycle
,enerd
,shakev_vir
,groups
,
1666 &state
->hist
,vsite
,fp_field
,ed
,force_flags
);
1670 /* Construct the virtual sites for the initial configuration */
1671 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,NULL
,
1672 top
->idef
.iparams
,top
->idef
.il
,
1673 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
1679 /* note reversed shakev_vir and shake_vir, this affects which ones are gathered in the step */
1680 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
1681 wcycle
,enerd
,force_vir
,shake_vir
,shakev_vir
,total_vir
,ekin
,pres
,mu_tot
,
1682 constr
,&chkpt
,&terminate
,&terminate_now
,&(nlh
.nabnsb
),state
->box
,
1683 top_global
,&pcurr
,&bSumEkinhOld
,bRerunMD
,bEkinFullStep
,FALSE
,
1684 TRUE
,FALSE
,FALSE
,FALSE
,FALSE
,TRUE
,Flags
& MD_READ_EKIN
,top_global
->natoms
);
1686 /* Calculate the initial half step temperature */
1687 temp0
= enerd
->term
[F_TEMP
];
1689 /* Initiate data for the special cases */
1692 bufstate
= init_bufstate(state
->natoms
,(ir
->opts
.ngtc
+1)*NNHCHAIN
); /* extra state for barostat */
1697 snew(xcopy
,state
->natoms
);
1698 snew(vcopy
,state
->natoms
);
1699 copy_rvecn(state
->x
,xcopy
,0,state
->natoms
);
1700 copy_rvecn(state
->v
,vcopy
,0,state
->natoms
);
1701 copy_mat(state
->box
,boxcopy
);
1706 /* need to make an initial call to get the Trotter Mass variables set. */
1707 trotter_update(ir
,ekind
,enerd
,state
,ekin
,total_vir
,mdatoms
,&MassQ
,FALSE
,FALSE
,FALSE
,FALSE
);
1712 if (constr
&& !ir
->bContinuation
&& ir
->eConstrAlg
== econtLINCS
)
1715 "RMS relative constraint deviation after constraining: %.2e\n",
1716 constr_rmsd(constr
,FALSE
));
1718 fprintf(fplog
,"Initial temperature: %g K\n",enerd
->term
[F_TEMP
]);
1721 fprintf(stderr
,"starting md rerun '%s', reading coordinates from"
1722 " input trajectory '%s'\n\n",
1723 *(top_global
->name
),opt2fn("-rerun",nfile
,fnm
));
1726 fprintf(stderr
,"Calculated time to finish depends on nsteps from "
1727 "run input file,\nwhich may not correspond to the time "
1728 "needed to process input trajectory.\n\n");
1733 fprintf(stderr
,"starting mdrun '%s'\n",
1734 *(top_global
->name
));
1735 if (ir
->init_step
> 0)
1737 fprintf(stderr
,"%s steps, %8.1f ps (continuing from step %s, %8.1f ps).\n",
1738 gmx_step_str(ir
->init_step
+ir
->nsteps
,sbuf
),
1739 (ir
->init_step
+ir
->nsteps
)*ir
->delta_t
,
1740 gmx_step_str(ir
->init_step
,sbuf2
),
1741 ir
->init_step
*ir
->delta_t
);
1745 fprintf(stderr
,"%s steps, %8.1f ps.\n",
1746 gmx_step_str(ir
->nsteps
,sbuf
),ir
->nsteps
*ir
->delta_t
);
1749 fprintf(fplog
,"\n");
1752 /* Set and write start time */
1753 runtime_start(runtime
);
1754 print_date_and_time(fplog
,cr
->nodeid
,"Started mdrun",runtime
);
1755 wallcycle_start(wcycle
,ewcRUN
);
1757 fprintf(fplog
,"\n");
1759 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
1761 chkpt_ret
=fcCheckPointParallel( cr
->nodeid
,
1763 if ( chkpt_ret
== 0 )
1764 gmx_fatal( 3,__FILE__
,__LINE__
, "Checkpoint error on step %d\n", 0 );
1768 /***********************************************************
1770 * Loop over MD steps
1772 ************************************************************/
1774 /* if rerunMD then read coordinates and velocities from input trajectory */
1777 if (getenv("GMX_FORCE_UPDATE"))
1779 bForceUpdate
= TRUE
;
1782 bNotLastFrame
= read_first_frame(oenv
,&status
,
1783 opt2fn("-rerun",nfile
,fnm
),
1784 &rerun_fr
,TRX_NEED_X
| TRX_READ_V
);
1785 if (rerun_fr
.natoms
!= top_global
->natoms
)
1788 "Number of atoms in trajectory (%d) does not match the "
1789 "run input file (%d)\n",
1790 rerun_fr
.natoms
,top_global
->natoms
);
1792 if (ir
->ePBC
!= epbcNONE
)
1796 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
);
1798 if (max_cutoff2(ir
->ePBC
,rerun_fr
.box
) < sqr(fr
->rlistlong
))
1800 gmx_fatal(FARGS
,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr
.step
,rerun_fr
.time
);
1803 /* Set the shift vectors.
1804 * Necessary here when have a static box different from the tpr box.
1806 calc_shifts(rerun_fr
.box
,fr
->shift_vec
);
1810 /* loop over MD steps or if rerunMD to end of input trajectory */
1812 /* Skip the first Nose-Hoover integration when we get the state from tpx */
1813 bStateFromTPX
= !opt2bSet("-cpi",nfile
,fnm
);
1814 bInitStep
= bFirstStep
&& (bStateFromTPX
|| ir
->eI
==eiVV
);
1816 bSumEkinhOld
= FALSE
;
1819 step
= ir
->init_step
;
1822 if (ir
->nstlist
== -1)
1824 init_nlistheuristics(&nlh
,bGStatEveryStep
,step
);
1827 bLastStep
= (bRerunMD
|| step_rel
> ir
->nsteps
);
1828 while (!bLastStep
|| (bRerunMD
&& bNotLastFrame
)) {
1830 wallcycle_start(wcycle
,ewcSTEP
);
1832 GMX_MPE_LOG(ev_timestep1
);
1835 if (rerun_fr
.bStep
) {
1836 step
= rerun_fr
.step
;
1837 step_rel
= step
- ir
->init_step
;
1839 if (rerun_fr
.bTime
) {
1849 bLastStep
= (step_rel
== ir
->nsteps
);
1850 t
= t0
+ step
*ir
->delta_t
;
1853 if (ir
->efep
!= efepNO
)
1855 if (bRerunMD
&& rerun_fr
.bLambda
&& (ir
->delta_lambda
!=0))
1857 state_global
->lambda
= rerun_fr
.lambda
;
1861 state_global
->lambda
= lam0
+ step
*ir
->delta_lambda
;
1863 state
->lambda
= state_global
->lambda
;
1864 bDoDHDL
= do_per_step(step
,ir
->nstdhdl
);
1869 update_annealing_target_temp(&(ir
->opts
),t
);
1874 if (!(DOMAINDECOMP(cr
) && !MASTER(cr
)))
1876 for(i
=0; i
<state_global
->natoms
; i
++)
1878 copy_rvec(rerun_fr
.x
[i
],state_global
->x
[i
]);
1882 for(i
=0; i
<state_global
->natoms
; i
++)
1884 copy_rvec(rerun_fr
.v
[i
],state_global
->v
[i
]);
1889 for(i
=0; i
<state_global
->natoms
; i
++)
1891 clear_rvec(state_global
->v
[i
]);
1895 fprintf(stderr
,"\nWARNING: Some frames do not contain velocities.\n"
1896 " Ekin, temperature and pressure are incorrect,\n"
1897 " the virial will be incorrect when constraints are present.\n"
1899 bRerunWarnNoV
= FALSE
;
1903 copy_mat(rerun_fr
.box
,state_global
->box
);
1904 copy_mat(state_global
->box
,state
->box
);
1906 if (vsite
&& (Flags
& MD_RERUN_VSITE
))
1908 if (DOMAINDECOMP(cr
))
1910 gmx_fatal(FARGS
,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
1914 /* Following is necessary because the graph may get out of sync
1915 * with the coordinates if we only have every N'th coordinate set
1917 mk_mshift(fplog
,graph
,fr
->ePBC
,state
->box
,state
->x
);
1918 shift_self(graph
,state
->box
,state
->x
);
1920 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,state
->v
,
1921 top
->idef
.iparams
,top
->idef
.il
,
1922 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
1925 unshift_self(graph
,state
->box
,state
->x
);
1930 /* Stop Center of Mass motion */
1931 bStopCM
= (ir
->comm_mode
!= ecmNO
&& do_per_step(step
,ir
->nstcomm
));
1933 /* Copy back starting coordinates in case we're doing a forcefield scan */
1936 for(ii
=0; (ii
<state
->natoms
); ii
++)
1938 copy_rvec(xcopy
[ii
],state
->x
[ii
]);
1939 copy_rvec(vcopy
[ii
],state
->v
[ii
]);
1941 copy_mat(boxcopy
,state
->box
);
1946 /* for rerun MD always do Neighbour Searching */
1947 bNS
= (bFirstStep
|| ir
->nstlist
!= 0);
1952 /* Determine whether or not to do Neighbour Searching and LR */
1953 bNStList
= (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0);
1955 bNS
= (bFirstStep
|| bExchanged
|| bNStList
||
1956 (ir
->nstlist
== -1 && nlh
.nabnsb
> 0));
1958 if (bNS
&& ir
->nstlist
== -1)
1960 set_nlistheuristics(&nlh
,bFirstStep
|| bExchanged
,step
);
1964 if (terminate_now
> 0 || (terminate_now
< 0 && bNS
))
1969 /* Determine whether or not to update the Born radii if doing GB */
1970 bBornRadii
=bFirstStep
;
1971 if (ir
->implicit_solvent
&& (step
% ir
->nstgbradii
==0))
1976 do_log
= do_per_step(step
,ir
->nstlog
) || bFirstStep
|| bLastStep
;
1977 do_verbose
= bVerbose
&&
1978 (step
% stepout
== 0 || bFirstStep
|| bLastStep
);
1980 if (bNS
&& !(bFirstStep
&& ir
->bContinuation
&& !bRerunMD
))
1984 bMasterState
= TRUE
;
1988 bMasterState
= FALSE
;
1989 /* Correct the new box if it is too skewed */
1990 if (DYNAMIC_BOX(*ir
))
1992 if (correct_box(fplog
,step
,state
->box
,graph
))
1994 bMasterState
= TRUE
;
1997 if (DOMAINDECOMP(cr
) && bMasterState
)
1999 dd_collect_state(cr
->dd
,state
,state_global
);
2003 if (DOMAINDECOMP(cr
))
2005 /* Repartition the domain decomposition */
2006 wallcycle_start(wcycle
,ewcDOMDEC
);
2007 dd_partition_system(fplog
,step
,cr
,
2008 bMasterState
,nstglobalcomm
,
2009 state_global
,top_global
,ir
,
2010 state
,&f
,mdatoms
,top
,fr
,
2011 vsite
,shellfc
,constr
,
2012 nrnb
,wcycle
,do_verbose
);
2013 wallcycle_stop(wcycle
,ewcDOMDEC
);
2017 if (MASTER(cr
) && do_log
&& !bFFscan
)
2019 print_ebin_header(fplog
,step
,t
,state
->lambda
);
2022 if (ir
->efep
!= efepNO
)
2024 update_mdatoms(mdatoms
,state
->lambda
);
2027 if (bRerunMD
&& rerun_fr
.bV
)
2030 /* We need the kinetic energy at minus the half step for determining
2031 * the full step kinetic energy and possibly for T-coupling.*/
2032 /* This is almost certainly not quite working correctly yet */
2033 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2034 wcycle
,enerd
,NULL
,NULL
,NULL
,NULL
,ekin
,NULL
,mu_tot
,
2035 constr
,&chkpt
,&terminate
,&terminate_now
,&(nlh
.nabnsb
),state
->box
,
2036 top_global
,&pcurr
,&bSumEkinhOld
,TRUE
,FALSE
,FALSE
,
2037 TRUE
,FALSE
,FALSE
,FALSE
,FALSE
,FALSE
,FALSE
,top_global
->natoms
);
2039 clear_mat(force_vir
);
2041 /* Ionize the atoms if necessary */
2044 ionize(fplog
,oenv
,mdatoms
,top_global
,t
,ir
,state
->x
,state
->v
,
2045 mdatoms
->start
,mdatoms
->start
+mdatoms
->homenr
,state
->box
,cr
);
2048 /* Update force field in ffscan program */
2051 if (update_forcefield(fplog
,
2053 mdatoms
->nr
,state
->x
,state
->box
)) {
2054 if (gmx_parallel_env())
2062 GMX_MPE_LOG(ev_timestep2
);
2064 if ((bNS
|| bLastStep
) && (step
> ir
->init_step
) && !bRerunMD
)
2066 bCPT
= (chkpt
> 0 || (bLastStep
&& (Flags
& MD_CONFOUT
)));
2077 /* Determine the pressure:
2078 * always when we want exact averages in the energy file,
2079 * at ns steps when we have pressure coupling,
2080 * otherwise only at energy output steps (set below).
2082 bNstEner
= (bGStatEveryStep
|| do_per_step(step
,ir
->nstcalcenergy
));
2083 bCalcEner
= bNstEner
;
2084 bCalcPres
= (bGStatEveryStep
|| (ir
->epc
!= epcNO
&& bNS
));
2086 /* Do we need global communication ? */
2087 bGStat
= (bCalcEner
|| bStopCM
||
2088 (ir
->nstlist
== -1 && !bRerunMD
&& step
>= nlh
.step_nscheck
));
2090 do_ene
= (do_per_step(step
,ir
->nstenergy
) || bLastStep
);
2092 if (do_ene
|| do_log
)
2099 force_flags
= (GMX_FORCE_STATECHANGED
|
2100 GMX_FORCE_ALLFORCES
|
2101 (bNStList
? GMX_FORCE_DOLR
: 0) |
2103 (bCalcEner
? GMX_FORCE_VIRIAL
: 0) |
2104 (bDoDHDL
? GMX_FORCE_DHDL
: 0));
2108 /* Now is the time to relax the shells */
2109 count
=relax_shell_flexcon(fplog
,cr
,bVerbose
,bFFscan
? step
+1 : step
,
2111 bStopCM
,top
,top_global
,
2113 state
,f
,force_vir
,mdatoms
,
2114 nrnb
,wcycle
,graph
,groups
,
2115 shellfc
,fr
,bBornRadii
,t
,mu_tot
,
2116 state
->natoms
,&bConverged
,vsite
,
2127 /* The coordinates (x) are shifted (to get whole molecules)
2129 * This is parallellized as well, and does communication too.
2130 * Check comments in sim_util.c
2133 do_force(fplog
,cr
,ir
,step
,nrnb
,wcycle
,top
,top_global
,groups
,
2134 state
->box
,state
->x
,&state
->hist
,
2135 f
,force_vir
,mdatoms
,enerd
,fcd
,
2136 state
->lambda
,graph
,
2137 fr
,vsite
,mu_tot
,t
,fp_field
,ed
,bBornRadii
,
2138 (bNS
? GMX_FORCE_NS
: 0) | force_flags
);
2141 GMX_BARRIER(cr
->mpi_comm_mygroup
);
2145 mu_aver
= calc_mu_aver(cr
,state
->x
,mdatoms
->chargeA
,
2146 mu_tot
,&top_global
->mols
,mdatoms
,gnx
,grpindex
);
2149 if (bTCR
&& bFirstStep
)
2151 tcr
=init_coupling(fplog
,nfile
,fnm
,cr
,fr
,mdatoms
,&(top
->idef
));
2152 fprintf(fplog
,"Done init_coupling\n");
2156 /* in the first case, we only have to iterate with Trotter-based pressure control */
2158 /* ############### START FIRST UPDATE HALF-STEP ############### */
2161 if ( !bRerunMD
|| bForceUpdate
) {
2162 // if ( !bRerunMD || rerun_fr.bV || bForceUpdate) { // Why was rerun_fr.bV here? Doesn't make sense.
2163 wallcycle_start(wcycle
,ewcUPDATE
);
2166 update_coords(fplog
,step
,ir
,mdatoms
,state
,
2167 f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2168 ekind
,NULL
,M
,wcycle
,upd
,bInitStep
,TRUE
);
2171 bIterate
= ((ir
->epc
== epcTROTTER
) && (constr
) && (!bInitStep
));
2172 bFirstIterate
= TRUE
;
2173 /* for iterations, we save these vectors, as we will be self-consistently iterating
2175 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
2177 /* save the state */
2179 copy_coupling_state(state
,bufstate
,NULL
,NULL
,
2180 ekin
,ekin_save
,ekind
,ekind_save
);
2183 while (bIterate
|| bFirstIterate
)
2187 copy_coupling_state(bufstate
,state
,NULL
,NULL
,
2188 ekin_save
,ekin
,ekind_save
,ekind
);
2192 if (bFirstIterate
&& bTrotter
)
2194 /* The first time, we need a decent first estimate
2195 of veta(t+dt) to compute the constraints. Do
2196 this by computing the box volume part of the
2197 trotter integration at this time. Nothing else
2198 should be changed by this routine here. If
2199 !(first time), we start with the previous value
2201 veta_save
= state
->veta
;
2202 trotter_update(ir
,ekind
,enerd
,state
,ekin
,total_vir
,mdatoms
,&MassQ
,TRUE
,FALSE
,TRUE
,bInitStep
);
2203 vetanew
= state
->veta
;
2204 state
->veta
= veta_save
;
2208 if ( !bRerunMD
|| bForceUpdate
)
2209 //if ( !bRerunMD || rerun_fr.bV || bForceUpdate) // Why was rerun_fr.bV here? Doesn't make sense.
2211 update_constraints(fplog
,step
,&dvdl
,ir
,mdatoms
,state
,graph
,f
,
2212 &top
->idef
,shakev_vir
,NULL
,
2213 cr
,nrnb
,wcycle
,upd
,constr
,
2214 bInitStep
,TRUE
,bCalcPres
,vetanew
);
2215 if (!bOK
&& !bFFscan
)
2217 gmx_fatal(FARGS
,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2219 if (fr
->bSepDVDL
&& fplog
&& do_log
)
2221 fprintf(fplog
,sepdvdlformat
,"Constraint",0.0,dvdl
);
2223 enerd
->term
[F_DHDL_CON
] += dvdl
;
2227 { /* Need to unshift here if a do_force has been
2228 called in the previous step */
2229 unshift_self(graph
,state
->box
,state
->x
);
2232 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2233 wcycle
,enerd
,force_vir
,shakev_vir
,shake_vir
,total_vir
,ekin
,pres
,mu_tot
,
2234 constr
,&chkpt
,&terminate
,&terminate_now
,&(nlh
.nabnsb
),state
->box
,
2235 (bTCR
&& bFFscan
)?top_global
:NULL
,&pcurr
,&bSumEkinhOld
,
2236 bRerunMD
,bEkinFullStep
,bDoCM
,bGStat
,bNEMD
,TRUE
,bIterate
,
2237 bFirstIterate
,FALSE
,FALSE
,top_global
->natoms
);
2239 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
2242 trotter_update(ir
,ekind
,enerd
,state
,ekin
,total_vir
,mdatoms
,&MassQ
,TRUE
,TRUE
,TRUE
,bInitStep
);
2245 /* iterate until veta is converged. This -should- be equivalent
2246 to self consistently solving for the Lagrange multipliers */
2249 fprintf(debug
,"Iterating vir, ekin: %20.12f %20.12f\n",trace(total_vir
),trace(ekin
));
2251 if (done_iterating(&bFirstIterate
,&bIterate
,state
->veta
,&vetanew
,0)) break;
2254 if (bTrotter
&& !bInitStep
)
2256 /* update temperature and kinetic energy now that step is over */
2257 enerd
->term
[F_TEMP
] = sum_ekin(FALSE
,&(ir
->opts
),ekind
,ekin
,NULL
,bEkinFullStep
);
2258 enerd
->term
[F_EKIN
] = trace(ekin
);
2261 wallcycle_stop(wcycle
,ewcUPDATE
);
2262 GMX_MPE_LOG(ev_timestep1
);
2264 /* MRS -- done iterating -- compute the conserved quantity */
2266 if (ir
->etc
== etcTROTTER
|| ir
->etc
== etcTROTTEREKINH
)
2269 NPT_energy(ir
,state
->nosehoover_xi
,state
->nosehoover_vxi
,
2270 state
->boxv
, state
->veta
, state
->box
, &MassQ
);
2271 last_conserved
+= enerd
->term
[F_EKIN
];
2272 if ((ir
->eDispCorr
!= edispcEnerPres
) && (ir
->eDispCorr
!= edispcAllEnerPres
))
2274 last_conserved
-= enerd
->term
[F_DISPCORR
];
2278 /* ######## END FIRST UPDATE STEP ############## */
2279 /* ######## If doing VV, we now have v(dt) ###### */
2281 /* ################## START TRAJECTORY OUTPUT ################# */
2283 /* Now we have the energies and forces corresponding to the
2284 * coordinates at time t. We must output all of this before
2286 * for RerunMD t is read from input trajectory
2288 GMX_MPE_LOG(ev_output_start
);
2290 bX
= do_per_step(step
,ir
->nstxout
);
2291 bV
= do_per_step(step
,ir
->nstvout
);
2292 bF
= do_per_step(step
,ir
->nstfout
);
2293 bXTC
= do_per_step(step
,ir
->nstxtcout
);
2297 fcReportProgress( ir
->nsteps
, step
);
2299 bX
= bX
|| bLastStep
; /*enforce writing positions and velocities
2301 bV
= bV
|| bLastStep
;
2303 int nthreads
=(cr
->nthreads
==0 ? 1 : cr
->nthreads
);
2304 int nnodes
=(cr
->nnodes
==0 ? 1 : cr
->nnodes
);
2307 /*Gromacs drives checkpointing; no ||
2308 fcCheckPointPendingThreads(cr->nodeid,
2310 /* sync bCPT and fc record-keeping */
2311 if (bCPT
&& MASTER(cr
))
2312 fcRequestCheckPoint();
2317 if (bX
|| bV
|| bF
|| bXTC
|| bCPT
)
2319 wallcycle_start(wcycle
,ewcTRAJ
);
2322 if (state
->flags
& (1<<estLD_RNG
))
2324 get_stochd_state(upd
,state
);
2330 state_global
->ekinstate
.bUpToDate
= FALSE
;
2334 update_ekinstate(&state_global
->ekinstate
,ekind
);
2335 state_global
->ekinstate
.bUpToDate
= TRUE
;
2337 update_energyhistory(&state_global
->enerhist
,mdebin
);
2340 write_traj(fplog
,cr
,fp_trn
,bX
,bV
,bF
,fp_xtc
,bXTC
,ir
->xtcprec
,
2341 fn_cpt
,bCPT
,top_global
,ir
->eI
,ir
->simulation_part
,
2342 step
,t
,state
,state_global
,f
,f_global
,&n_xtc
,&x_xtc
);
2344 if (bLastStep
&& step_rel
== ir
->nsteps
&&
2345 (Flags
& MD_CONFOUT
) && MASTER(cr
) &&
2346 !bRerunMD
&& !bFFscan
)
2348 /* x and v have been collected in write_traj */
2349 fprintf(stderr
,"\nWriting final coordinates.\n");
2350 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
&&
2353 /* Make molecules whole only for confout writing */
2354 do_pbc_mtop(fplog
,ir
->ePBC
,state
->box
,top_global
,state_global
->x
);
2356 write_sto_conf_mtop(ftp2fn(efSTO
,nfile
,fnm
),
2357 *top_global
->name
,top_global
,
2358 state_global
->x
,state_global
->v
,
2359 ir
->ePBC
,state
->box
);
2362 wallcycle_stop(wcycle
,ewcTRAJ
);
2364 GMX_MPE_LOG(ev_output_finish
);
2366 /* ################## END TRAJECTORY OUTPUT ################ */
2368 /* Determine the pressure:
2369 * always when we want exact averages in the energy file,
2370 * at ns steps when we have pressure coupling,
2371 * otherwise only at energy output steps (set below).
2374 bNstEner
= (bGStatEveryStep
|| do_per_step(step
,ir
->nstcalcenergy
));
2375 bCalcEner
= bNstEner
;
2376 bCalcPres
= (bGStatEveryStep
|| (ir
->epc
!= epcNO
&& bNS
));
2378 /* Do we need global communication ? */
2379 bGStat
= (bGStatEveryStep
|| bStopCM
|| bNS
||
2380 (ir
->nstlist
== -1 && !bRerunMD
&& step
>= nlh
.step_nscheck
));
2382 do_ene
= (do_per_step(step
,ir
->nstenergy
) || bLastStep
);
2384 if (do_ene
|| do_log
)
2390 bIterate
= ((ir
->epc
== epcTROTTER
) && (constr
));
2391 bFirstIterate
= TRUE
;
2393 /* for iterations, we save these vectors, as we will be redoing the calculations */
2396 copy_coupling_state(state
,bufstate
,NULL
,NULL
,ekin
,ekin_save
,ekind
,ekind_save
);
2399 while (bIterate
|| bFirstIterate
)
2401 /* We now restore these vectors to redo the calculation with improved extended variables */
2404 copy_coupling_state(bufstate
,state
,NULL
,NULL
,ekin_save
,ekin
,ekind_save
,ekind
);
2407 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
2408 so scroll down for that logic */
2410 /* ######### START SECOND UPDATE STEP ################# */
2412 GMX_MPE_LOG(ev_update_start
);
2414 if (!bRerunMD
|| rerun_fr
.bV
|| bForceUpdate
)
2416 wallcycle_start(wcycle
,ewcUPDATE
);
2420 /* Box is changed in update() when we do pressure coupling,
2421 * but we should still use the old box for energy corrections and when
2422 * writing it to the energy file, so it matches the trajectory files for
2423 * the same timestep above. Make a copy in a separate array.
2426 copy_mat(state
->box
,lastbox
);
2429 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS*/
2434 m_add(shakev_vir
,force_vir
,total_vir
);
2435 if (bFirstIterate
&& bInitStep
)
2437 /* Double shakevir to approximate shake_vir for a better initial guess the first time.
2438 for all subsequent searches, we just start with the last virial. Should be a better
2439 way to do this . .. */
2440 m_add(shakev_vir
,total_vir
,total_vir
);
2444 if (abs(trace(shake_vir
))!=0)
2446 scalevir
= tracevir
/trace(shake_vir
);
2452 msmul(shake_vir
,scalevir
,shake_vir
);
2453 m_add(total_vir
,shake_vir
,total_vir
);
2456 trotter_update(ir
,ekind
,enerd
,state
,ekin
,total_vir
,mdatoms
,&MassQ
,FALSE
,TRUE
,TRUE
,bInitStep
);
2459 /* We can only do Berendsen coupling after we have summed the kinetic
2460 * energy or virial. Since the happens in global_state after update,
2461 * we should only do it at step % nstlist = 1 with bGStatEveryStep=FALSE.
2464 update_extended(fplog
,step
,ir
,mdatoms
,state
,ekind
,ekin
,pcoupl_mu
,M
,wcycle
,
2465 upd
,bInitStep
,FALSE
,&MassQ
);
2467 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2468 ekind
,ekin
,M
,wcycle
,upd
,bInitStep
,FALSE
);
2470 update_constraints(fplog
,step
,&dvdl
,ir
,mdatoms
,state
,graph
,f
,
2471 &top
->idef
,shake_vir
,force_vir
,
2472 cr
,nrnb
,wcycle
,upd
,constr
,
2473 bInitStep
,FALSE
,bCalcPres
,state
->veta
); /* !!!examine this!!! */
2475 if (!bOK
&& !bFFscan
)
2477 gmx_fatal(FARGS
,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2480 if (fr
->bSepDVDL
&& fplog
&& do_log
)
2482 fprintf(fplog
,sepdvdlformat
,"Constraint",0.0,dvdl
);
2484 enerd
->term
[F_DHDL_CON
] += dvdl
;
2485 wallcycle_stop(wcycle
,ewcUPDATE
);
2489 /* Need to unshift here */
2490 unshift_self(graph
,state
->box
,state
->x
);
2493 GMX_BARRIER(cr
->mpi_comm_mygroup
);
2494 GMX_MPE_LOG(ev_update_finish
);
2498 wallcycle_start(wcycle
,ewcVSITECONSTR
);
2501 shift_self(graph
,state
->box
,state
->x
);
2503 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,state
->v
,
2504 top
->idef
.iparams
,top
->idef
.il
,
2505 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
2509 unshift_self(graph
,state
->box
,state
->x
);
2511 wallcycle_stop(wcycle
,ewcVSITECONSTR
);
2514 /* ############## IF NOT VV, CALCULATE EKIN AND PRESSURE HERE ############ */
2515 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2516 wcycle
,enerd
,force_vir
,shake_vir
,shakev_vir
,total_vir
,ekin
,pres
,mu_tot
,
2517 constr
,&chkpt
,&terminate
,&terminate_now
,&(nlh
.nabnsb
),lastbox
,
2518 (bTCR
&& bFFscan
)?top_global
:NULL
,&pcurr
,&bSumEkinhOld
,
2519 bRerunMD
,bEkinFullStep
,bDoCM
,bGStat
,bNEMD
,FALSE
,bIterate
,
2520 bFirstIterate
,FALSE
,FALSE
,top_global
->natoms
);
2522 /* ############# END CALC EKIN AND PRESSURE ################# */
2524 if (done_iterating(&bFirstIterate
,&bIterate
,trace(shake_vir
),&tracevir
,1)) break;
2527 update_box(fplog
,step
,ir
,mdatoms
,state
,graph
,f
,
2528 scale_tot
,pcoupl_mu
,nrnb
,wcycle
,upd
,bInitStep
,FALSE
);
2531 /* ################# END UPDATE STEP 2 ################# */
2532 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
2534 /* Determine the wallclock run time up till now */
2535 run_time
= (double)time(NULL
) - (double)runtime
->real
;
2537 /* Check whether everything is still allright */
2538 if ((bGotTermSignal
|| bGotUsr1Signal
) && !bHandledSignal
)
2540 if (bGotTermSignal
|| ir
->nstlist
== 0)
2550 terminate_now
= terminate
;
2555 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2556 bGotTermSignal
? "TERM" : "USR1",
2557 terminate
==-1 ? "NS " : "");
2561 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2562 bGotTermSignal
? "TERM" : "USR1",
2563 terminate
==-1 ? "NS " : "");
2565 bHandledSignal
=TRUE
;
2567 else if (MASTER(cr
) && (bNS
|| ir
->nstlist
<= 0) &&
2568 (max_hours
> 0 && run_time
> max_hours
*60.0*60.0*0.99) &&
2571 /* Signal to terminate the run */
2572 terminate
= (ir
->nstlist
== 0 ? 1 : -1);
2575 terminate_now
= terminate
;
2579 fprintf(fplog
,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step
,sbuf
),max_hours
*0.99);
2581 fprintf(stderr
, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step
,sbuf
),max_hours
*0.99);
2584 if (ir
->nstlist
== -1 && !bRerunMD
)
2586 /* When bGStatEveryStep=FALSE, global_stat is only called
2587 * when we check the atom displacements, not at NS steps.
2588 * This means that also the bonded interaction count check is not
2589 * performed immediately after NS. Therefore a few MD steps could
2590 * be performed with missing interactions.
2591 * But wrong energies are never written to file,
2592 * since energies are only written after global_stat
2595 if (step
>= nlh
.step_nscheck
)
2597 nlh
.nabnsb
= natoms_beyond_ns_buffer(ir
,fr
,&top
->cgs
,
2598 nlh
.scale_tot
,state
->x
);
2602 /* This is not necessarily true,
2603 * but step_nscheck is determined quite conservatively.
2609 /* In parallel we only have to check for checkpointing in steps
2610 * where we do global communication,
2611 * otherwise the other nodes don't know.
2613 if (MASTER(cr
) && ((bGStat
|| !PAR(cr
)) &&
2616 run_time
>= nchkpt
*cpt_period
*60.0)))
2625 /* The coordinates (x) were unshifted in update */
2626 if (bFFscan
&& (shellfc
==NULL
|| bConverged
))
2628 if (print_forcefield(fplog
,enerd
->term
,mdatoms
->homenr
,
2630 &(top_global
->mols
),mdatoms
->massT
,pres
))
2632 if (gmx_parallel_env()) {
2635 fprintf(stderr
,"\n");
2642 /* Only do GCT when the relaxation of shells (minimization) has converged,
2643 * otherwise we might be coupling to bogus energies.
2644 * In parallel we must always do this, because the other sims might
2648 /* Since this is called with the new coordinates state->x, I assume
2649 * we want the new box state->box too. / EL 20040121
2651 do_coupling(fplog
,oenv
,nfile
,fnm
,tcr
,t
,step
,enerd
->term
,fr
,
2653 mdatoms
,&(top
->idef
),mu_aver
,
2654 top_global
->mols
.nr
,cr
,
2655 state
->box
,total_vir
,pres
,
2656 mu_tot
,state
->x
,f
,bConverged
);
2660 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
2662 sum_dhdl(enerd
,state
->lambda
,ir
);
2663 enerd
->term
[F_ETOT
] = enerd
->term
[F_EPOT
] + enerd
->term
[F_EKIN
];
2672 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_EPOT
] + last_conserved
;
2674 case etcTROTTEREKINH
:
2675 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] + last_conserved
;
2678 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] +
2679 NPT_energy(ir
,state
->nosehoover_xi
,
2680 state
->nosehoover_vxi
,state
->boxv
,state
->veta
,state
->box
,&MassQ
);
2683 enerd
->term
[F_ECONSERVED
] =
2684 /* not the best name to be using here -- think about renaming later */
2685 enerd
->term
[F_ETOT
] + vrescale_energy(&(ir
->opts
),
2686 state
->therm_integral
);
2692 /* Check for excessively large energies */
2696 real etot_max
= 1e200
;
2698 real etot_max
= 1e30
;
2700 if (fabs(enerd
->term
[F_ETOT
]) > etot_max
)
2702 fprintf(stderr
,"Energy too large (%g), giving up\n",
2703 enerd
->term
[F_ETOT
]);
2706 /* ######### END PREPARING EDR OUTPUT ########### */
2708 /* Time for performance */
2709 if (((step
% stepout
) == 0) || bLastStep
)
2711 runtime_upd_proc(runtime
);
2721 upd_mdebin(mdebin
,bDoDHDL
? fp_dhdl
: NULL
,TRUE
,
2722 t
,mdatoms
->tmass
,enerd
,state
,lastbox
,
2723 shake_vir
,force_vir
,total_vir
,pres
,
2724 ekind
,mu_tot
,constr
);
2728 upd_mdebin_step(mdebin
);
2731 do_dr
= do_per_step(step
,ir
->nstdisreout
);
2732 do_or
= do_per_step(step
,ir
->nstorireout
);
2734 print_ebin(fp_ene
,do_ene
,do_dr
,do_or
,do_log
?fplog
:NULL
,step
,t
,
2735 eprNORMAL
,bCompact
,mdebin
,fcd
,groups
,&(ir
->opts
));
2737 if (ir
->ePull
!= epullNO
)
2739 pull_print_output(ir
->pull
,step
,t
);
2742 if (do_per_step(step
,ir
->nstlog
))
2744 if(fflush(fplog
) != 0)
2746 gmx_fatal(FARGS
,"Cannot flush logfile - maybe you are out of quota?");
2752 /* Remaining runtime */
2753 if (MULTIMASTER(cr
) && do_verbose
)
2757 fprintf(stderr
,"\n");
2759 print_time(stderr
,runtime
,step
,ir
);
2762 /* Replica exchange */
2764 if ((repl_ex_nst
> 0) && (step
> 0) && !bLastStep
&&
2765 do_per_step(step
,repl_ex_nst
))
2767 bExchanged
= replica_exchange(fplog
,cr
,repl_ex
,
2768 state_global
,enerd
->term
,
2771 if (bExchanged
&& PAR(cr
))
2773 if (DOMAINDECOMP(cr
))
2775 dd_partition_system(fplog
,step
,cr
,TRUE
,1,
2776 state_global
,top_global
,ir
,
2777 state
,&f
,mdatoms
,top
,fr
,
2778 vsite
,shellfc
,constr
,
2783 bcast_state(cr
,state
,FALSE
);
2790 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
2791 /* Complicated conditional when bGStatEveryStep=FALSE.
2792 * We can not just use bGStat, since then the simulation results
2793 * would depend on nstenergy and nstlog or step_nscheck.
2795 if ((state
->flags
& (1<<estPRES_PREV
)) &&
2797 (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0) ||
2798 (ir
->nstlist
< 0 && nlh
.nabnsb
> 0) ||
2799 (ir
->nstlist
== 0 && bGStat
)))
2801 /* Store the pressure in t_state for pressure coupling
2802 * at the next MD step.
2804 copy_mat(pres
,state
->pres_prev
);
2807 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
2811 /* read next frame from input trajectory */
2812 bNotLastFrame
= read_next_frame(oenv
,status
,&rerun_fr
);
2815 if (!bRerunMD
|| !rerun_fr
.bStep
)
2817 /* increase the MD step number */
2822 cycles
= wallcycle_stop(wcycle
,ewcSTEP
);
2823 if (DOMAINDECOMP(cr
) && wcycle
)
2825 dd_cycles_add(cr
->dd
,cycles
,ddCyclStep
);
2828 if (step_rel
== wcycle_get_reset_counters(wcycle
))
2830 /* Reset all the counters related to performance over the run */
2831 reset_all_counters(fplog
,cr
,step
,&step_rel
,ir
,wcycle
,nrnb
,runtime
);
2832 wcycle_set_reset_counters(wcycle
, 0);
2835 /* End of main MD loop */
2839 runtime_end(runtime
);
2846 if (!(cr
->duty
& DUTY_PME
))
2848 /* Tell the PME only node to finish */
2854 if (bGStatEveryStep
&& !bRerunMD
)
2856 print_ebin(fp_ene
,FALSE
,FALSE
,FALSE
,fplog
,step
,t
,
2857 eprAVER
,FALSE
,mdebin
,fcd
,groups
,&(ir
->opts
));
2867 gmx_fio_fclose(fp_dhdl
);
2871 gmx_fio_fclose(fp_field
);
2876 if (ir
->nstlist
== -1 && nlh
.nns
> 0 && fplog
)
2878 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
)));
2879 fprintf(fplog
,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh
.ab
/nlh
.nns
);
2882 if (shellfc
&& fplog
)
2884 fprintf(fplog
,"Fraction of iterations that converged: %.2f %%\n",
2885 (nconverged
*100.0)/step_rel
);
2886 fprintf(fplog
,"Average number of force evaluations per MD step: %.2f\n\n",
2890 if (repl_ex_nst
> 0 && MASTER(cr
))
2892 print_replica_exchange_statistics(fplog
,repl_ex
);
2895 runtime
->nsteps_done
= step_rel
;