4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Copyright (c) 1991-2001, University of Groningen, The Netherlands
12 * This program is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU General Public License
14 * as published by the Free Software Foundation; either version 2
15 * of the License, or (at your option) any later version.
17 * If you want to redistribute modifications, please consider that
18 * scientific software is very special. Version control is crucial -
19 * bugs must be traceable. We will be happy to consider code for
20 * inclusion in the official distribution, but derived work must not
21 * be called official GROMACS. Details are found in the README & COPYING
22 * files - if they are missing, get the official version at www.gromacs.org.
24 * To help us fund GROMACS development, we humbly ask that you cite
25 * the papers on the package - you can find them in the top README file.
27 * For more info, check our website at http://www.gromacs.org
30 * GROningen Mixture of Alchemy and Childrens' Stories
32 static char *SRCID_md_c
= "$Id$";
65 volatile bool bGotTermSignal
= FALSE
, bGotUsr1Signal
= FALSE
;
67 static RETSIGTYPE
signal_handler(int n
)
71 bGotTermSignal
= TRUE
;
74 bGotUsr1Signal
= TRUE
;
81 void mdrunner(t_commrec
*cr
,t_commrec
*mcr
,int nfile
,t_filenm fnm
[],
82 bool bVerbose
,bool bCompact
,
83 int nDlb
,int nstepout
,t_edsamyn
*edyn
,
86 double nodetime
=0,realtime
;
88 rvec
*buf
,*f
,*vold
,*v
,*vt
,*x
,box_size
;
100 bool bDummies
,bParDummies
;
101 t_comm_dummies dummycomm
;
105 /* Initiate everything (snew sets to zero!) */
112 snew(nrnb
,cr
->nnodes
);
114 if (bVerbose
&& MASTER(cr
))
115 fprintf(stderr
,"Getting Loaded...\n");
118 /* The master thread on the master node reads from disk, then passes everything
119 * around the ring, and finally frees the stuff
122 distribute_parts(cr
->left
,cr
->right
,cr
->nodeid
,cr
->nnodes
,parm
,
123 ftp2fn(efTPX
,nfile
,fnm
),nDlb
);
125 /* Every node (including the master) reads the data from the ring */
126 init_parts(stdlog
,cr
,
127 parm
,top
,&x
,&v
,&mdatoms
,nsb
,
128 MASTER(cr
) ? LIST_SCALARS
| LIST_PARM
: 0,
129 &bParDummies
,&dummycomm
);
132 init_single(stdlog
,parm
,ftp2fn(efTPX
,nfile
,fnm
),top
,&x
,&v
,&mdatoms
,nsb
);
135 snew(buf
,nsb
->natoms
);
137 snew(vt
,nsb
->natoms
);
138 snew(vold
,nsb
->natoms
);
140 if (bVerbose
&& MASTER(cr
))
141 fprintf(stderr
,"Loaded with Money\n\n");
143 /* Index numbers for parallellism... */
144 nsb
->nodeid
= cr
->nodeid
;
145 top
->idef
.nodeid
= cr
->nodeid
;
147 /* Group stuff (energies etc) */
148 init_groups(stdlog
,mdatoms
,&(parm
->ir
.opts
),grps
);
149 /* Copy the cos acceleration to the groups struct */
150 grps
->cosacc
.cos_accel
= parm
->ir
.cos_accel
;
152 /* Periodicity stuff */
153 graph
=mk_graph(&(top
->idef
),top
->atoms
.nr
,FALSE
,FALSE
);
155 p_graph(debug
,"Initial graph",graph
);
157 /* Distance Restraints */
158 init_disres(stdlog
,top
->idef
.il
[F_DISRES
].nr
,top
->idef
.il
[F_DISRES
].iatoms
,
159 top
->idef
.iparams
,&(parm
->ir
),mcr
,fcd
);
161 /* Orientation restraints */
162 init_orires(stdlog
,top
->idef
.il
[F_ORIRES
].nr
,top
->idef
.il
[F_ORIRES
].iatoms
,
163 top
->idef
.iparams
,x
,mdatoms
,&(parm
->ir
),mcr
,fcd
);
165 /* check if there are dummies */
167 for(i
=0; (i
<F_NRE
) && !bDummies
; i
++)
168 bDummies
= ((interaction_function
[i
].flags
& IF_DUMMY
) &&
169 (top
->idef
.il
[i
].nr
> 0));
171 /* Initiate forcerecord */
173 init_forcerec(stdlog
,fr
,&(parm
->ir
),top
,cr
,mdatoms
,nsb
,parm
->box
,FALSE
,
174 opt2fn("-table",nfile
,fnm
),FALSE
);
175 fr
->bSepDVDL
= ((Flags
& MD_SEPDVDL
) == MD_SEPDVDL
);
177 for(m
=0; (m
<DIM
); m
++)
178 box_size
[m
]=parm
->box
[m
][m
];
180 /* Initiate PPPM if necessary */
181 if (fr
->eeltype
== eelPPPM
)
182 init_pppm(stdlog
,cr
,nsb
,FALSE
,TRUE
,box_size
,getenv("GMXGHAT"),&parm
->ir
);
183 if (fr
->eeltype
== eelPME
)
184 init_pme(stdlog
,cr
,parm
->ir
.nkx
,parm
->ir
.nky
,parm
->ir
.nkz
,parm
->ir
.pme_order
,
185 HOMENR(nsb
),parm
->ir
.bOptFFT
,parm
->ir
.ewald_geometry
);
187 /* Now do whatever the user wants us to do (how flexible...) */
188 switch (parm
->ir
.eI
) {
192 start_t
=do_md(stdlog
,cr
,mcr
,nfile
,fnm
,
193 bVerbose
,bCompact
,bDummies
,
194 bParDummies
? &dummycomm
: NULL
,
195 nstepout
,parm
,grps
,top
,ener
,fcd
,x
,vold
,v
,vt
,f
,buf
,
196 mdatoms
,nsb
,nrnb
,graph
,edyn
,fr
,box_size
,Flags
);
199 start_t
=do_cg(stdlog
,nfile
,fnm
,parm
,top
,grps
,nsb
,
200 x
,f
,buf
,mdatoms
,parm
->ekin
,ener
,fcd
,
201 nrnb
,bVerbose
,bDummies
,
202 bParDummies
? &dummycomm
: NULL
,
203 cr
,mcr
,graph
,fr
,box_size
);
206 start_t
=do_steep(stdlog
,nfile
,fnm
,parm
,top
,grps
,nsb
,
207 x
,f
,buf
,mdatoms
,parm
->ekin
,ener
,fcd
,
208 nrnb
,bVerbose
,bDummies
,
209 bParDummies
? &dummycomm
: NULL
,
210 cr
,mcr
,graph
,fr
,box_size
);
213 start_t
=do_nm(stdlog
,cr
,nfile
,fnm
,
214 bVerbose
,bCompact
,nstepout
,parm
,grps
,
215 top
,ener
,fcd
,x
,vold
,v
,vt
,f
,buf
,
216 mdatoms
,nsb
,nrnb
,graph
,edyn
,fr
,box_size
);
219 fatal_error(0,"Invalid integrator (%d)...\n",parm
->ir
.eI
);
222 /* Some timing stats */
224 realtime
=difftime(time(NULL
),start_t
);
225 if ((nodetime
=node_time()) == 0)
231 /* Convert back the atoms */
232 md2atoms(mdatoms
,&(top
->atoms
),TRUE
);
234 /* Finish up, write some stuff
235 * if rerunMD, don't write last frame again
237 finish_run(stdlog
,cr
,ftp2fn(efSTO
,nfile
,fnm
),
238 nsb
,top
,parm
,nrnb
,nodetime
,realtime
,parm
->ir
.nsteps
,
239 parm
->ir
.eI
==eiMD
|| parm
->ir
.eI
==eiSD
|| parm
->ir
.eI
==eiBD
);
241 /* Does what it says */
242 print_date_and_time(stdlog
,cr
->nodeid
,"Finished mdrun");
245 time_t do_md(FILE *log
,t_commrec
*cr
,t_commrec
*mcr
,int nfile
,t_filenm fnm
[],
246 bool bVerbose
,bool bCompact
,
247 bool bDummies
, t_comm_dummies
*dummycomm
,
248 int stepout
,t_parm
*parm
,t_groups
*grps
,t_topology
*top
,
249 real ener
[],t_fcdata
*fcd
,
250 rvec x
[],rvec vold
[],rvec v
[],rvec vt
[],rvec f
[],
251 rvec buf
[],t_mdatoms
*mdatoms
,t_nsborder
*nsb
,t_nrnb nrnb
[],
252 t_graph
*graph
,t_edsamyn
*edyn
,t_forcerec
*fr
,rvec box_size
,
256 int fp_ene
=0,fp_trn
=0,step
;
259 real t
,lambda
,t0
,lam0
,SAfactor
;
260 bool bNS
,bStopCM
,bTYZ
,bRerunMD
,bNotLastFrame
=FALSE
,
261 bFirstStep
,bLastStep
,bNEMD
,do_log
,bRerunWarnNoV
=TRUE
;
262 tensor force_vir
,pme_vir
,shake_vir
;
264 char *traj
,*xtc_traj
; /* normal and compressed trajectory filename */
270 t_pull pulldata
; /* for pull code */
271 /* A boolean (disguised as a real) to terminate mdrun */
274 /* XMDRUN stuff: shell, general coupling etc. */
276 int nshell
,nshell_tot
,count
,nconverged
=0;
277 t_shell
*shells
=NULL
;
280 bool bShell
,bIonize
=FALSE
,bGlas
=FALSE
;
281 bool bTCR
=FALSE
,bConverged
=TRUE
,bOK
;
286 t_coupl_rec
*tcr
=NULL
;
288 /* End of XMDRUN stuff */
290 /* Turn on signal handling */
291 signal(SIGTERM
,signal_handler
);
292 signal(SIGUSR1
,signal_handler
);
294 /* Check for special mdrun options */
295 bRerunMD
= (Flags
& MD_RERUN
) == MD_RERUN
;
296 bIonize
= (Flags
& MD_IONIZE
) == MD_IONIZE
;
297 bGlas
= (Flags
& MD_GLAS
) == MD_GLAS
;
298 bFFscan
= (Flags
& MD_FFSCAN
) == MD_FFSCAN
;
301 init_md(cr
,&parm
->ir
,parm
->box
,&t
,&t0
,&lambda
,&lam0
,&SAfactor
,
303 nfile
,fnm
,&traj
,&xtc_traj
,&fp_ene
,&fp_dgdl
,&mdebin
,grps
,
304 force_vir
,pme_vir
,shake_vir
,mdatoms
,mu_tot
,&bNEMD
,&vcm
,nsb
);
307 /* Check for polarizable models */
308 shells
= init_shells(log
,START(nsb
),HOMENR(nsb
),&top
->idef
,
312 gmx_sumi(1,&nshell_tot
,cr
);
313 bShell
= nshell_tot
> 0;
315 /* Check whether we have to do dipole stuff */
316 if (opt2bSet("-dn",nfile
,fnm
))
317 rd_index(opt2fn("-dn",nfile
,fnm
),1,&gnx
,&grpindex
,&grpname
);
319 gnx
= top
->blocks
[ebMOLS
].nr
;
321 for(i
=0; (i
<gnx
); i
++)
325 /* Check whether we have to GCT stuff */
326 bTCR
= ftp2bSet(efGCT
,nfile
,fnm
);
327 if (MASTER(cr
) && bTCR
)
328 fprintf(stderr
,"Will do General Coupling Theory!\n");
330 /* Remove periodicity */
331 if (fr
->ePBC
!= epbcNONE
)
332 do_pbc_first(log
,parm
,box_size
,fr
,graph
,x
);
335 /* Initialize pull code */
336 init_pull(log
,nfile
,fnm
,&pulldata
,x
,mdatoms
,parm
->box
,START(nsb
), HOMENR(nsb
), cr
);
337 if (pulldata
.bPull
&& PAR(cr
) && pulldata
.runtype
==eConstraint
)
338 fatal_error(0,"Can not do constraint runs in parallel");
340 if (!parm
->ir
.bUncStart
)
341 do_shakefirst(log
,bTYZ
,lambda
,ener
,parm
,nsb
,mdatoms
,x
,vold
,buf
,f
,v
,
342 graph
,cr
,&mynrnb
,grps
,fr
,top
,edyn
,&pulldata
);
345 /* Compute initial EKin for all.. */
346 if (grps
->cosacc
.cos_accel
== 0)
347 calc_ke_part(TRUE
,parm
->ir
.eI
==eiSD
,
348 START(nsb
),HOMENR(nsb
),vold
,v
,vt
,&(parm
->ir
.opts
),
349 mdatoms
,grps
,&mynrnb
,lambda
,&ener
[F_DVDLKIN
]);
351 calc_ke_part_visc(TRUE
,START(nsb
),HOMENR(nsb
),
352 parm
->box
,x
,vold
,v
,vt
,&(parm
->ir
.opts
),
353 mdatoms
,grps
,&mynrnb
,lambda
,&ener
[F_DVDLKIN
]);
357 global_stat(log
,cr
,ener
,force_vir
,shake_vir
,
358 &(parm
->ir
.opts
),grps
,&mynrnb
,nrnb
,vcm
,&terminate
);
361 /* Calculate Temperature coupling parameters lambda */
362 ener
[F_TEMP
] = sum_ekin(&(parm
->ir
.opts
),grps
,parm
->ekin
,bTYZ
);
363 if(parm
->ir
.etc
==etcBERENDSEN
)
364 berendsen_tcoupl(&(parm
->ir
.opts
),grps
,
365 parm
->ir
.delta_t
,SAfactor
);
366 else if(parm
->ir
.etc
==etcNOSEHOOVER
)
367 nosehoover_tcoupl(&(parm
->ir
.opts
),grps
,
368 parm
->ir
.delta_t
,SAfactor
);
371 /* Initiate data for the special cases */
373 snew(xcopy
,nsb
->natoms
);
374 for(ii
=0; (ii
<nsb
->natoms
); ii
++)
375 copy_rvec(x
[ii
],xcopy
[ii
]);
378 /* Write start time and temperature */
379 start_t
=print_date_and_time(log
,cr
->nodeid
,"Started mdrun");
382 fprintf(log
,"Initial temperature: %g K\n",ener
[F_TEMP
]);
384 fprintf(stderr
,"starting md rerun '%s', reading coordinates from"
385 " input trajectory '%s'\n\n",
386 *(top
->name
),opt2fn("-rerun",nfile
,fnm
));
388 fprintf(stderr
,"Calculated time to finish depends on nsteps from "
389 "run input file,\nwhich may not correspond to the time "
390 "needed to process input trajectory.\n\n");
392 fprintf(stderr
,"starting mdrun '%s'\n%d steps, %8.1f ps.\n\n",
393 *(top
->name
),parm
->ir
.nsteps
,parm
->ir
.nsteps
*parm
->ir
.delta_t
);
395 /* Set the node time counter to 0 after initialisation */
398 /***********************************************************
402 ************************************************************/
404 /* if rerunMD then read coordinates and velocities from input trajectory */
406 bNotLastFrame
= read_first_frame(&status
,opt2fn("-rerun",nfile
,fnm
),
407 &rerun_fr
,TRX_NEED_X
| TRX_READ_V
);
408 if (rerun_fr
.natoms
!= mdatoms
->nr
)
409 fatal_error(0,"Number of atoms in trajectory (%d) does not match the "
410 "run input file (%d)\n",rerun_fr
.natoms
,mdatoms
->nr
);
413 /* loop over MD steps or if rerunMD to end of input trajectory */
416 while ((!bRerunMD
&& (step
<=parm
->ir
.nsteps
)) ||
417 (bRerunMD
&& bNotLastFrame
)) {
421 step
= rerun_fr
.step
;
426 for(i
=0; i
<mdatoms
->nr
; i
++)
427 copy_rvec(rerun_fr
.x
[i
],x
[i
]);
429 for(i
=0; i
<mdatoms
->nr
; i
++)
430 copy_rvec(rerun_fr
.v
[i
],v
[i
]);
432 for(i
=0; i
<mdatoms
->nr
; i
++)
435 fprintf(stderr
,"\nWARNING: Some frames do not contain velocities.\n"
436 " Ekin, temperature and pressure are incorrect,\n"
437 " the virial will be incorrect when constraints are present.\n"
439 bRerunWarnNoV
= FALSE
;
442 copy_mat(rerun_fr
.box
,parm
->box
);
444 /* for rerun MD always do Neighbour Searching */
445 bNS
= ((parm
->ir
.nstlist
!=0) || bFirstStep
);
447 /* Determine whether or not to do Neighbour Searching */
448 bNS
= ((parm
->ir
.nstlist
&& (step
% parm
->ir
.nstlist
==0)) || bFirstStep
);
450 bLastStep
=(step
==parm
->ir
.nsteps
);
452 do_log
= do_per_step(step
,parm
->ir
.nstlog
) || bLastStep
;
454 /* Stop Center of Mass motion */
455 bStopCM
= do_per_step(step
,abs(parm
->ir
.nstcomm
));
457 /* Copy back starting coordinates in case we're doing a forcefield scan */
459 for(ii
=0; (ii
<nsb
->natoms
); ii
++)
460 copy_rvec(xcopy
[ii
],x
[ii
]);
464 shift_self(graph
,parm
->box
,x
);
466 construct_dummies(log
,x
,&mynrnb
,parm
->ir
.delta_t
,v
,&top
->idef
,
467 graph
,cr
,parm
->box
,dummycomm
);
469 unshift_self(graph
,parm
->box
,x
);
474 /* Set values for invmass etc. This routine not parallellized, but hardly
475 * ever used, only when doing free energy calculations.
477 init_mdatoms(mdatoms
,lambda
,bFirstStep
);
479 clear_mat(force_vir
);
481 /* Ionize the atoms if necessary */
483 ionize(log
,mdatoms
,top
->atoms
.atomname
,t
,&parm
->ir
,x
,v
,
484 START(nsb
),START(nsb
)+HOMENR(nsb
),parm
->box
,cr
);
486 /* Update force field in ffscan program */
488 update_forcefield(nfile
,fnm
,fr
);
491 /* Now is the time to relax the shells */
492 count
=relax_shells(log
,cr
,mcr
,bVerbose
,step
,parm
,bNS
,bStopCM
,top
,
493 ener
,fcd
,x
,vold
,v
,vt
,f
,buf
,mdatoms
,nsb
,&mynrnb
,graph
,
494 grps
,force_vir
,pme_vir
,bShell
,
495 nshell
,shells
,fr
,traj
,t
,lambda
,mu_tot
,
496 nsb
->natoms
,parm
->box
,&bConverged
);
503 /* The coordinates (x) are shifted (to get whole molecules) in do_force
504 * This is parallellized as well, and does communication too.
505 * Check comments in sim_util.c
507 do_force(log
,cr
,mcr
,parm
,nsb
,force_vir
,pme_vir
,step
,&mynrnb
,
508 top
,grps
,x
,v
,f
,buf
,mdatoms
,ener
,fcd
,bVerbose
&& !PAR(cr
),
509 lambda
,graph
,bNS
,FALSE
,fr
,mu_tot
,FALSE
);
513 mu_aver
=calc_mu_aver(cr
,nsb
,x
,mdatoms
->chargeA
,mu_tot
,top
,mdatoms
,
517 do_glas(log
,START(nsb
),HOMENR(nsb
),x
,f
,
518 fr
,mdatoms
,top
->idef
.atnr
,&parm
->ir
,ener
);
520 if (bTCR
&& (step
== 0)) {
521 tcr
=init_coupling(log
,nfile
,fnm
,cr
,fr
,mdatoms
,&(top
->idef
));
522 fprintf(log
,"Done init_coupling\n");
526 /* Now we have the energies and forces corresponding to the
527 * coordinates at time t. We must output all of this before
529 * for RerunMD t is read from input trajectory
532 t
= t0
+ step
*parm
->ir
.delta_t
;
534 if (parm
->ir
.efep
!= efepNO
) {
535 if (bRerunMD
&& rerun_fr
.bLambda
)
536 lambda
= rerun_fr
.lambda
;
538 lambda
= lam0
+ step
*parm
->ir
.delta_lambda
;
540 if (parm
->ir
.bSimAnn
) {
541 SAfactor
= 1.0 - t
/parm
->ir
.zero_temp_time
;
546 if (MASTER(cr
) && do_log
&& !bFFscan
)
547 print_ebin_header(log
,step
,t
,lambda
,SAfactor
);
550 spread_dummy_f(log
,x
,f
,&mynrnb
,&top
->idef
,dummycomm
,cr
);
552 /* Calculation of the virial must be done after dummies! */
553 /* Question: Is it correct to do the PME forces after this? */
554 calc_virial(log
,START(nsb
),HOMENR(nsb
),x
,f
,
555 force_vir
,pme_vir
,graph
,parm
->box
,&mynrnb
,fr
,FALSE
);
557 /* Spread the LR force on dummy particle to the other particles...
558 * This is parallellized. MPI communication is performed
559 * if the constructing atoms aren't local.
561 if (bDummies
&& fr
->bEwald
)
562 spread_dummy_f(log
,x
,fr
->f_pme
,&mynrnb
,&top
->idef
,dummycomm
,cr
);
564 sum_lrforces(f
,fr
,START(nsb
),HOMENR(nsb
));
566 xx
= (do_per_step(step
,parm
->ir
.nstxout
) || bLastStep
) ? x
: NULL
;
567 vv
= (do_per_step(step
,parm
->ir
.nstvout
) || bLastStep
) ? v
: NULL
;
568 ff
= (do_per_step(step
,parm
->ir
.nstfout
)) ? f
: NULL
;
570 fp_trn
= write_traj(log
,cr
,traj
,nsb
,step
,t
,lambda
,
571 nrnb
,nsb
->natoms
,xx
,vv
,ff
,parm
->box
);
574 /* don't write xtc and last structure for rerunMD */
575 if (!bRerunMD
&& !bFFscan
) {
576 if (do_per_step(step
,parm
->ir
.nstxtcout
))
577 write_xtc_traj(log
,cr
,xtc_traj
,nsb
,mdatoms
,
578 step
,t
,x
,parm
->box
,parm
->ir
.xtcprec
);
579 if (bLastStep
&& MASTER(cr
)) {
580 fprintf(stderr
,"Writing final coordinates.\n");
581 write_sto_conf(ftp2fn(efSTO
,nfile
,fnm
),
582 *top
->name
, &(top
->atoms
),x
,v
,parm
->box
);
586 clear_mat(shake_vir
);
588 /* All pulling except constraint type pulling happens before the update,
589 * Constraint happens in update
591 if(pulldata
.bPull
&& !(pulldata
.runtype
== eConstraint
))
592 pull(&pulldata
,x
,f
,parm
->box
,top
,parm
->ir
.delta_t
,step
,
593 mdatoms
->nr
,mdatoms
, START(nsb
), HOMENR(nsb
), cr
);
596 clear_rvecs(nsb
->natoms
,buf
);
598 /* This is also parallellized, but check code in update.c */
599 /* bOK = update(nsb->natoms,START(nsb),HOMENR(nsb),step,lambda,&ener[F_DVDL], */
601 update(nsb
->natoms
,START(nsb
),HOMENR(nsb
),step
,lambda
,&ener
[F_DVDL
],
602 parm
,SAfactor
,mdatoms
,x
,graph
,f
,buf
,vold
,vt
,v
,
603 top
,grps
,shake_vir
,cr
,&mynrnb
,bTYZ
,TRUE
,edyn
,&pulldata
,bNEMD
);
604 if (!bOK
&& !bFFscan
)
605 fatal_error(0,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
607 /* The coordinates (x) were unshifted in update */
608 if (bFFscan
&& (!bShell
|| bConverged
))
609 print_forcefield(log
,ener
[F_EPOT
],HOMENR(nsb
),f
,buf
,xcopy
,
610 &(top
->blocks
[ebMOLS
]),mdatoms
->massT
);
612 if (parm
->ir
.epc
!=epcNO
)
613 correct_box(parm
->box
,fr
,graph
);
614 /* (un)shifting should NOT be done after this,
615 * since the box vectors might have changed
618 /* Non-equilibrium MD: this is parallellized, but only does communication
619 * when there really is NEMD.
621 if (PAR(cr
) && bNEMD
)
622 accumulate_u(cr
,&(parm
->ir
.opts
),grps
);
625 if (grps
->cosacc
.cos_accel
== 0)
626 calc_ke_part(FALSE
,parm
->ir
.eI
==eiSD
,
627 START(nsb
),HOMENR(nsb
),vold
,v
,vt
,&(parm
->ir
.opts
),
628 mdatoms
,grps
,&mynrnb
,lambda
,&ener
[F_DVDLKIN
]);
630 calc_ke_part_visc(FALSE
,START(nsb
),HOMENR(nsb
),
631 parm
->box
,x
,vold
,v
,vt
,&(parm
->ir
.opts
),
632 mdatoms
,grps
,&mynrnb
,lambda
,&ener
[F_DVDLKIN
]);
635 /* Calculate center of mass velocity if necessary, also parallellized */
636 if (bStopCM
&& !bFFscan
)
637 calc_vcm_grp(log
,START(nsb
),HOMENR(nsb
),mdatoms
->massT
,x
,v
,vcm
);
639 /* Check whether everything is still allright */
640 if (bGotTermSignal
|| bGotUsr1Signal
) {
645 fprintf(log
,"\n\nReceived the %s signal\n\n",
646 bGotTermSignal
? "TERM" : "USR1");
649 fprintf(stderr
,"\n\nReceived the %s signal\n\n",
650 bGotTermSignal
? "TERM" : "USR1");
653 bGotTermSignal
= FALSE
;
654 bGotUsr1Signal
= FALSE
;
658 /* Globally (over all NODEs) sum energy, virial etc.
659 * This includes communication
661 global_stat(log
,cr
,ener
,force_vir
,shake_vir
,
662 &(parm
->ir
.opts
),grps
,&mynrnb
,nrnb
,vcm
,&terminate
);
663 /* Correct for double counting energies, should be moved to
666 if (fr
->bTwinRange
&& !bNS
)
667 for(i
=0; (i
<grps
->estat
.nn
); i
++) {
668 grps
->estat
.ee
[egLR
][i
] /= cr
->nnodes
;
669 grps
->estat
.ee
[egLJLR
][i
] /= cr
->nnodes
;
673 cp_nrnb(&(nrnb
[0]),&mynrnb
);
675 /* This is just for testing. Nothing is actually done to Ekin
676 * since that would require extra communication.
678 if (!bNEMD
&& debug
&& (vcm
->nr
> 0))
679 correct_ekin(debug
,START(nsb
),START(nsb
)+HOMENR(nsb
),v
,
681 mdatoms
->massT
,mdatoms
->tmass
,parm
->ekin
);
683 if ((terminate
!= 0) && (step
< parm
->ir
.nsteps
)) {
684 if (terminate
<0 && parm
->ir
.nstxout
)
685 /* this is the USR1 signal and we are writing x to trr,
686 stop at next x frame in trr */
687 parm
->ir
.nsteps
= (step
/ parm
->ir
.nstxout
+ 1) * parm
->ir
.nstxout
;
689 parm
->ir
.nsteps
= step
+1;
690 fprintf(log
,"\nSetting nsteps to %d\n\n",parm
->ir
.nsteps
);
693 fprintf(stderr
,"\nSetting nsteps to %d\n\n",parm
->ir
.nsteps
);
696 /* erase the terminate signal */
700 /* Do center of mass motion removal */
701 if (bStopCM
&& !bFFscan
) {
702 check_cm_grp(log
,vcm
);
703 do_stopcm_grp(log
,START(nsb
),HOMENR(nsb
),x
,v
,vcm
);
704 inc_nrnb(&mynrnb
,eNR_STOPCM
,HOMENR(nsb
));
706 calc_vcm_grp(log,START(nsb),HOMENR(nsb),mdatoms->massT,x,v,vcm);
707 check_cm_grp(log,vcm);
708 do_stopcm_grp(log,START(nsb),HOMENR(nsb),x,v,vcm);
709 check_cm_grp(log,vcm);
713 /* Add force and shake contribution to the virial */
714 m_add(force_vir
,shake_vir
,parm
->vir
);
716 /* Sum the potential energy terms from group contributions */
717 sum_epot(&(parm
->ir
.opts
),grps
,ener
);
719 /* Calculate the amplitude of the cosine velocity profile */
720 grps
->cosacc
.vcos
= grps
->cosacc
.mvcos
/mdatoms
->tmass
;
722 /* Sum the kinetic energies of the groups & calc temp */
723 ener
[F_TEMP
]=sum_ekin(&(parm
->ir
.opts
),grps
,parm
->ekin
,bTYZ
);
724 ener
[F_EKIN
]=trace(parm
->ekin
);
725 ener
[F_ETOT
]=ener
[F_EPOT
]+ener
[F_EKIN
];
727 /* Check for excessively large energies */
728 if (bIonize
&& fabs(ener
[F_ETOT
]) > 1e18
) {
729 fprintf(stderr
,"Energy too large (%g), giving up\n",ener
[F_ETOT
]);
733 /* Calculate Temperature coupling parameters lambda and adjust
734 * target temp when doing simulated annealing
736 if(parm
->ir
.etc
==etcBERENDSEN
)
737 berendsen_tcoupl(&(parm
->ir
.opts
),grps
,parm
->ir
.delta_t
,SAfactor
);
738 else if(parm
->ir
.etc
==etcNOSEHOOVER
)
739 nosehoover_tcoupl(&(parm
->ir
.opts
),grps
,parm
->ir
.delta_t
,SAfactor
);
741 /* Calculate pressure and apply LR correction if PPPM is used */
742 calc_pres(fr
->ePBC
,parm
->box
,parm
->ekin
,parm
->vir
,parm
->pres
,
743 (fr
->eeltype
==eelPPPM
) ? ener
[F_LR
] : 0.0);
745 /* Calculate long range corrections to pressure and energy */
747 set_avcsix(log
,fr
,mdatoms
);
749 /* Calculate long range corrections to pressure and energy */
750 calc_dispcorr(log
,parm
->ir
.eDispCorr
,
751 fr
,mdatoms
->nr
,parm
->box
,parm
->pres
,parm
->vir
,ener
);
754 /* Only do GCT when the relaxation of shells (minimization) has converged,
755 * otherwise we might be coupling to bogus energies.
756 * In parallel we must always do this, because the other sims might
759 do_coupling(log
,nfile
,fnm
,tcr
,t
,step
,ener
,fr
,
760 &(parm
->ir
),MASTER(cr
),
761 mdatoms
,&(top
->idef
),mu_aver
,
762 top
->blocks
[ebMOLS
].nr
,(mcr
!=NULL
) ? mcr
: cr
,
763 parm
->box
,parm
->vir
,parm
->pres
,
764 mu_tot
,x
,f
,bConverged
);
768 /* Time for performance */
769 if (((step
% 10) == 0) || bLastStep
)
773 if (MASTER(cr
) && !bFFscan
) {
774 bool do_ene
,do_dr
,do_or
;
776 upd_mdebin(mdebin
,fp_dgdl
,mdatoms
->tmass
,step
,t
,ener
,parm
->box
,shake_vir
,
777 force_vir
,parm
->vir
,parm
->pres
,grps
,mu_tot
,
778 (parm
->ir
.etc
==etcNOSEHOOVER
));
779 do_ene
= do_per_step(step
,parm
->ir
.nstenergy
) || bLastStep
;
780 do_dr
= do_per_step(step
,parm
->ir
.nstdisreout
) || bLastStep
;
781 do_or
= do_per_step(step
,parm
->ir
.nstorireout
) || bLastStep
;
782 print_ebin(fp_ene
,do_ene
,do_dr
,do_or
,do_log
?log
:NULL
,step
,t
,
783 eprNORMAL
,bCompact
,mdebin
,fcd
,&(top
->atoms
));
788 /* Remaining runtime */
789 if (MASTER(cr
) && bVerbose
&& ( ((step
% stepout
)==0) || bLastStep
)) {
791 fprintf(stderr
,"\n");
792 print_time(stderr
,start_t
,step
,&parm
->ir
);
798 /* read next frame from input trajectory */
799 bNotLastFrame
= read_next_frame(status
,&rerun_fr
);
801 if (!bRerunMD
|| !rerun_fr
.bStep
)
802 /* increase the MD step number */
805 /* End of main MD loop */
808 /* Dump the NODE time to the log file on each node */
810 double *ct
,ctmax
,ctsum
;
813 ct
[cr
->nodeid
] = node_time();
814 gmx_sumd(cr
->nnodes
,ct
,cr
);
817 for(i
=1; (i
<cr
->nodeid
); i
++) {
818 ctmax
= max(ctmax
,ct
[i
]);
822 fprintf(log
,"\nTotal NODE time on node %d: %g\n",cr
->nodeid
,ct
[cr
->nodeid
]);
823 fprintf(log
,"Average NODE time: %g\n",ctsum
);
824 fprintf(log
,"Load imbalance reduced performance to %3d%% of max\n",
825 (int) (100.0*ctmax
/ctsum
));
829 print_ebin(fp_ene
,FALSE
,FALSE
,FALSE
,log
,step
,t
,
830 eprAVER
,FALSE
,mdebin
,fcd
,&(top
->atoms
));
831 print_ebin(fp_ene
,FALSE
,FALSE
,FALSE
,log
,step
,t
,
832 eprRMS
,FALSE
,mdebin
,fcd
,&(top
->atoms
));
834 if (!bRerunMD
&& parm
->ir
.nstxtcout
)
843 fprintf(log
,"Fraction of iterations that converged: %.2f %%\n",
844 (nconverged
*100.0)/(parm
->ir
.nsteps
+1));
845 fprintf(log
,"Average number of force evaluations per MD step: %.2f\n",
846 tcount
/(parm
->ir
.nsteps
+1));