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 * GROningen Mixture of Alchemy and Childrens' Stories
44 #include "gmx_fatal.h"
67 #include "checkpoint.h"
71 typedef struct gmx_global_stat
78 gmx_global_stat_t global_stat_init(t_inputrec *ir)
85 snew(gs->itc0,ir->opts.ngtc);
86 snew(gs->itc1,ir->opts.ngtc);
91 void global_stat_destroy(gmx_global_stat_t gs)
99 static int filter_enerdterm(real *afrom, gmx_bool bToBuffer, real *ato,
100 gmx_bool bTemp, gmx_bool bPres, gmx_bool bEner) {
105 for (i=0;i<F_NRE;i++)
121 ato[to++] = afrom[from++];
129 ato[to++] = afrom[from++];
135 ato[to++] = afrom[from++];
144 void global_stat(FILE *fplog,gmx_global_stat_t gs,
145 t_commrec *cr,gmx_enerdata_t *enerd,
146 tensor fvir,tensor svir,rvec mu_tot,
147 t_inputrec *inputrec,
148 gmx_ekindata_t *ekind,gmx_constr_t constr,
151 gmx_mtop_t *top_global, t_state *state_local,
152 gmx_bool bSumEkinhOld, int flags)
153 /* instead of current system, gmx_booleans for summing virial, kinetic energy, and other terms */
157 int ie=0,ifv=0,isv=0,irmsd=0,imu=0;
158 int idedl=0,idvdll=0,idvdlnl=0,iepl=0,icm=0,imass=0,ica=0,inb=0;
160 int icj=-1,ici=-1,icx=-1;
162 real copyenerd[F_NRE];
164 real *rmsd_data=NULL;
166 gmx_bool bVV,bTemp,bEner,bPres,bConstrVir,bEkinAveVel,bFirstIterate,bReadEkin;
168 bVV = EI_VV(inputrec->eI);
169 bTemp = flags & CGLO_TEMPERATURE;
170 bEner = flags & CGLO_ENERGY;
171 /* FIX ME after 4.5 */
172 /* temporary hack because we are using gmx_bool (unsigned char) */
173 bPres = (flags & CGLO_PRESSURE) != 0;
174 bConstrVir = (flags & CGLO_CONSTRAINT) != 0;
175 bFirstIterate = (flags & CGLO_FIRSTITERATE) != 0;
176 bEkinAveVel = (inputrec->eI==eiVV || (inputrec->eI==eiVVAK && bPres));
177 bReadEkin = (flags & CGLO_READEKIN) != 0;
185 /* This routine copies all the data to be summed to one big buffer
186 * using the t_bin struct.
189 /* First, we neeed to identify which enerd->term should be
190 communicated. Temperature and pressure terms should only be
191 communicated and summed when they need to be, to avoid repeating
192 the sums and overcounting. */
194 nener = filter_enerdterm(enerd->term,TRUE,copyenerd,bTemp,bPres,bEner);
196 /* First, the data that needs to be communicated with velocity verlet every time
197 This is just the constraint virial.*/
199 isv = add_binr(rb,DIM*DIM,svir[0]);
203 /* We need the force virial and the kinetic energy for the first time through with velocity verlet */
208 for(j=0; (j<inputrec->opts.ngtc); j++)
212 itc0[j]=add_binr(rb,DIM*DIM,ekind->tcstat[j].ekinh_old[0]);
214 if (bEkinAveVel && !bReadEkin)
216 itc1[j]=add_binr(rb,DIM*DIM,ekind->tcstat[j].ekinf[0]);
220 itc1[j]=add_binr(rb,DIM*DIM,ekind->tcstat[j].ekinh[0]);
223 /* these probably need to be put into one of these categories */
225 idedl = add_binr(rb,1,&(ekind->dekindl));
227 ica = add_binr(rb,1,&(ekind->cosacc.mvcos));
233 if ((bPres || !bVV) && bFirstIterate)
235 ifv = add_binr(rb,DIM*DIM,fvir[0]);
244 ie = add_binr(rb,nener,copyenerd);
249 rmsd_data = constr_rmsd_data(constr);
252 irmsd = add_binr(rb,inputrec->eI==eiSD2 ? 3 : 2,rmsd_data);
255 if (!NEED_MUTOT(*inputrec))
257 imu = add_binr(rb,DIM,mu_tot);
263 for(j=0; (j<egNR); j++)
265 inn[j]=add_binr(rb,enerd->grpp.nener,enerd->grpp.ener[j]);
268 if (inputrec->efep != efepNO)
270 idvdll = add_bind(rb,1,&enerd->dvdl_lin);
271 idvdlnl = add_bind(rb,1,&enerd->dvdl_nonlin);
272 if (enerd->n_lambda > 0)
274 iepl = add_bind(rb,enerd->n_lambda,enerd->enerpart_lambda);
281 icm = add_binr(rb,DIM*vcm->nr,vcm->group_p[0]);
283 imass = add_binr(rb,vcm->nr,vcm->group_mass);
285 if (vcm->mode == ecmANGULAR)
287 icj = add_binr(rb,DIM*vcm->nr,vcm->group_j[0]);
289 icx = add_binr(rb,DIM*vcm->nr,vcm->group_x[0]);
291 ici = add_binr(rb,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
296 if (DOMAINDECOMP(cr))
298 nb = cr->dd->nbonded_local;
299 inb = add_bind(rb,1,&nb);
304 isig = add_binr(rb,nsig,sig);
307 /* Global sum it all */
310 fprintf(debug,"Summing %d energies\n",rb->maxreal);
315 /* Extract all the data locally */
319 extract_binr(rb,isv ,DIM*DIM,svir[0]);
322 /* We need the force virial and the kinetic energy for the first time through with velocity verlet */
327 for(j=0; (j<inputrec->opts.ngtc); j++)
331 extract_binr(rb,itc0[j],DIM*DIM,ekind->tcstat[j].ekinh_old[0]);
333 if (bEkinAveVel && !bReadEkin) {
334 extract_binr(rb,itc1[j],DIM*DIM,ekind->tcstat[j].ekinf[0]);
338 extract_binr(rb,itc1[j],DIM*DIM,ekind->tcstat[j].ekinh[0]);
341 extract_binr(rb,idedl,1,&(ekind->dekindl));
342 extract_binr(rb,ica,1,&(ekind->cosacc.mvcos));
346 if ((bPres || !bVV) && bFirstIterate)
348 extract_binr(rb,ifv ,DIM*DIM,fvir[0]);
355 extract_binr(rb,ie,nener,copyenerd);
358 extract_binr(rb,irmsd,inputrec->eI==eiSD2 ? 3 : 2,rmsd_data);
360 if (!NEED_MUTOT(*inputrec))
362 extract_binr(rb,imu,DIM,mu_tot);
365 for(j=0; (j<egNR); j++)
367 extract_binr(rb,inn[j],enerd->grpp.nener,enerd->grpp.ener[j]);
369 if (inputrec->efep != efepNO)
371 extract_bind(rb,idvdll ,1,&enerd->dvdl_lin);
372 extract_bind(rb,idvdlnl,1,&enerd->dvdl_nonlin);
373 if (enerd->n_lambda > 0)
375 extract_bind(rb,iepl,enerd->n_lambda,enerd->enerpart_lambda);
378 /* should this be here, or with ekin?*/
381 extract_binr(rb,icm,DIM*vcm->nr,vcm->group_p[0]);
383 extract_binr(rb,imass,vcm->nr,vcm->group_mass);
385 if (vcm->mode == ecmANGULAR)
387 extract_binr(rb,icj,DIM*vcm->nr,vcm->group_j[0]);
389 extract_binr(rb,icx,DIM*vcm->nr,vcm->group_x[0]);
391 extract_binr(rb,ici,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
395 if (DOMAINDECOMP(cr))
397 extract_bind(rb,inb,1,&nb);
398 if ((int)(nb + 0.5) != cr->dd->nbonded_global)
400 dd_print_missing_interactions(fplog,cr,(int)(nb + 0.5),top_global,state_local);
405 filter_enerdterm(copyenerd,FALSE,enerd->term,bTemp,bPres,bEner);
406 /* Small hack for temp only - not entirely clear if still needed?*/
407 /* enerd->term[F_TEMP] /= (cr->nnodes - cr->npmenodes); */
413 extract_binr(rb,isig,nsig,sig);
418 int do_per_step(gmx_large_int_t step,gmx_large_int_t nstep)
421 return ((step % nstep)==0);
426 static void moveit(t_commrec *cr,
427 int left,int right,const char *s,rvec xx[])
432 move_rvecs(cr,FALSE,FALSE,left,right,
433 xx,NULL,(cr->nnodes-cr->npmenodes)-1,NULL);
436 gmx_mdoutf_t *init_mdoutf(int nfile,const t_filenm fnm[],int mdrun_flags,
437 const t_commrec *cr,const t_inputrec *ir,
438 const output_env_t oenv)
442 gmx_bool bAppendFiles;
452 of->eIntegrator = ir->eI;
453 of->simulation_part = ir->simulation_part;
457 bAppendFiles = (mdrun_flags & MD_APPENDFILES) != 0;
459 of->bKeepAndNumCPT = (mdrun_flags & MD_KEEPANDNUMCPT) != 0;
461 sprintf(filemode, bAppendFiles ? "a+" : "w+");
466 !(EI_DYNAMICS(ir->eI) &&
473 of->fp_trn = open_trn(ftp2fn(efTRN,nfile,fnm), filemode);
475 if (!EI_ENERGY_MINIMIZATION(ir->eI) &&
478 of->fp_xtc = open_xtc(ftp2fn(efXTC,nfile,fnm), filemode);
479 of->xtc_prec = ir->xtcprec;
481 of->fp_ene = open_enx(ftp2fn(efEDR,nfile,fnm), filemode);
482 of->fn_cpt = opt2fn("-cpo",nfile,fnm);
484 if (ir->efep != efepNO && ir->nstdhdl > 0 &&
485 (ir->separate_dhdl_file == sepdhdlfileYES ) &&
486 !EI_ENERGY_MINIMIZATION(ir->eI))
490 of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl",nfile,fnm),filemode);
494 of->fp_dhdl = open_dhdl(opt2fn("-dhdl",nfile,fnm),ir,oenv);
498 if (opt2bSet("-field",nfile,fnm) &&
499 (ir->ex[XX].n || ir->ex[YY].n || ir->ex[ZZ].n))
503 of->fp_dhdl = gmx_fio_fopen(opt2fn("-field",nfile,fnm),
508 of->fp_field = xvgropen(opt2fn("-field",nfile,fnm),
509 "Applied electric field","Time (ps)",
518 void done_mdoutf(gmx_mdoutf_t *of)
520 if (of->fp_ene != NULL)
522 close_enx(of->fp_ene);
526 close_xtc(of->fp_xtc);
530 close_trn(of->fp_trn);
532 if (of->fp_dhdl != NULL)
534 gmx_fio_fclose(of->fp_dhdl);
536 if (of->fp_field != NULL)
538 gmx_fio_fclose(of->fp_field);
544 void write_traj(FILE *fplog,t_commrec *cr,
547 gmx_mtop_t *top_global,
548 gmx_large_int_t step,double t,
549 t_state *state_local,t_state *state_global,
550 rvec *f_local,rvec *f_global,
551 int *n_xtc,rvec **x_xtc)
554 gmx_groups_t *groups;
559 #define MX(xvf) moveit(cr,GMX_LEFT,GMX_RIGHT,#xvf,xvf)
561 /* MRS -- defining these variables is to manage the difference
562 * between half step and full step velocities, but there must be a better way . . . */
564 local_v = state_local->v;
565 global_v = state_global->v;
567 if (DOMAINDECOMP(cr))
569 if (mdof_flags & MDOF_CPT)
571 dd_collect_state(cr->dd,state_local,state_global);
575 if (mdof_flags & (MDOF_X | MDOF_XTC))
577 dd_collect_vec(cr->dd,state_local,state_local->x,
580 if (mdof_flags & MDOF_V)
582 dd_collect_vec(cr->dd,state_local,local_v,
586 if (mdof_flags & MDOF_F)
588 dd_collect_vec(cr->dd,state_local,f_local,f_global);
593 if (mdof_flags & MDOF_CPT)
595 /* All pointers in state_local are equal to state_global,
596 * but we need to copy the non-pointer entries.
598 state_global->lambda = state_local->lambda;
599 state_global->veta = state_local->veta;
600 state_global->vol0 = state_local->vol0;
601 copy_mat(state_local->box,state_global->box);
602 copy_mat(state_local->boxv,state_global->boxv);
603 copy_mat(state_local->svir_prev,state_global->svir_prev);
604 copy_mat(state_local->fvir_prev,state_global->fvir_prev);
605 copy_mat(state_local->pres_prev,state_global->pres_prev);
609 /* Particle decomposition, collect the data on the master node */
610 if (mdof_flags & MDOF_CPT)
612 if (state_local->flags & (1<<estX)) MX(state_global->x);
613 if (state_local->flags & (1<<estV)) MX(state_global->v);
614 if (state_local->flags & (1<<estSDX)) MX(state_global->sd_X);
615 if (state_global->nrngi > 1) {
616 if (state_local->flags & (1<<estLD_RNG)) {
618 MPI_Gather(state_local->ld_rng ,
619 state_local->nrng*sizeof(state_local->ld_rng[0]),MPI_BYTE,
620 state_global->ld_rng,
621 state_local->nrng*sizeof(state_local->ld_rng[0]),MPI_BYTE,
622 MASTERRANK(cr),cr->mpi_comm_mygroup);
625 if (state_local->flags & (1<<estLD_RNGI))
628 MPI_Gather(state_local->ld_rngi,
629 sizeof(state_local->ld_rngi[0]),MPI_BYTE,
630 state_global->ld_rngi,
631 sizeof(state_local->ld_rngi[0]),MPI_BYTE,
632 MASTERRANK(cr),cr->mpi_comm_mygroup);
639 if (mdof_flags & (MDOF_X | MDOF_XTC)) MX(state_global->x);
640 if (mdof_flags & MDOF_V) MX(global_v);
642 if (mdof_flags & MDOF_F) MX(f_global);
648 if (mdof_flags & MDOF_CPT)
650 write_checkpoint(of->fn_cpt,of->bKeepAndNumCPT,
651 fplog,cr,of->eIntegrator,
652 of->simulation_part,step,t,state_global);
655 if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
657 fwrite_trn(of->fp_trn,step,t,state_local->lambda,
658 state_local->box,top_global->natoms,
659 (mdof_flags & MDOF_X) ? state_global->x : NULL,
660 (mdof_flags & MDOF_V) ? global_v : NULL,
661 (mdof_flags & MDOF_F) ? f_global : NULL);
662 if (gmx_fio_flush(of->fp_trn) != 0)
664 gmx_file("Cannot write trajectory; maybe you are out of quota?");
666 gmx_fio_check_file_position(of->fp_trn);
668 if (mdof_flags & MDOF_XTC) {
669 groups = &top_global->groups;
673 for(i=0; (i<top_global->natoms); i++)
675 if (ggrpnr(groups,egcXTC,i) == 0)
680 if (*n_xtc != top_global->natoms)
685 if (*n_xtc == top_global->natoms)
687 xxtc = state_global->x;
693 for(i=0; (i<top_global->natoms); i++)
695 if (ggrpnr(groups,egcXTC,i) == 0)
697 copy_rvec(state_global->x[i],xxtc[j++]);
701 if (write_xtc(of->fp_xtc,*n_xtc,step,t,
702 state_local->box,xxtc,of->xtc_prec) == 0)
704 gmx_fatal(FARGS,"XTC error - maybe you are out of quota?");
706 gmx_fio_check_file_position(of->fp_xtc);