3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
51 #include "checkpoint.h"
76 #include "mpelogging.h"
84 #include "checkpoint.h"
85 #include "mtop_util.h"
88 #include "sighandler.h"
101 /* We use the same defines as in mvdata.c here */
102 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d),(cr))
103 #define nblock_bc(cr,nr,d) gmx_bcast((nr)*sizeof((d)[0]), (d),(cr))
104 #define snew_bc(cr,d,nr) { if (!MASTER(cr)) snew((d),(nr)); }
106 /* The following two variables and the signal_handler function
107 * is used from pme.c as well
136 int natoms; /*nr of atoms per lipid*/
137 int mol1; /*id of the first lipid molecule*/
158 int search_string(char *s,int ng,char ***gn)
162 for(i=0; (i<ng); i++)
163 if (gmx_strcasecmp(s,*gn[i]) == 0)
166 gmx_fatal(FARGS,"Group %s not found in indexfile.\nMaybe you have non-default groups in your mdp file, while not using the '-n' option of grompp.\nIn that case use the '-n' option.\n",s);
171 int get_mol_id(int at,int nmblock,gmx_molblock_t *mblock, int *type, int *block)
176 for(i=0;i<nmblock;i++)
178 if(at<(mblock[i].nmol*mblock[i].natoms_mol))
180 mol_id+=at/mblock[i].natoms_mol;
181 *type = mblock[i].type;
185 at-= mblock[i].nmol*mblock[i].natoms_mol;
186 mol_id+=mblock[i].nmol;
190 gmx_fatal(FARGS,"Something is wrong in mol ids, at %d, mol_id %d",at,mol_id);
195 int get_block(int mol_id,int nmblock,gmx_molblock_t *mblock)
200 for(i=0;i<nmblock;i++)
202 nmol+=mblock[i].nmol;
207 gmx_fatal(FARGS,"mol_id %d larger than total number of molecules %d.\n",mol_id,nmol);
212 int get_tpr_version(const char *infile)
219 fio = open_tpx(infile,"r");
220 gmx_fio_checktype(fio);
222 precision = sizeof(real);
224 gmx_fio_do_string(fio,buf);
225 if (strncmp(buf,"VERSION",7))
226 gmx_fatal(FARGS,"Can not read file %s,\n"
227 " this file is from a Gromacs version which is older than 2.0\n"
228 " Make a new one with grompp or use a gro or pdb file, if possible",
229 gmx_fio_getname(fio));
230 gmx_fio_do_int(fio,precision);
231 bDouble = (precision == sizeof(double));
232 if ((precision != sizeof(float)) && !bDouble)
233 gmx_fatal(FARGS,"Unknown precision in file %s: real is %d bytes "
234 "instead of %d or %d",
235 gmx_fio_getname(fio),precision,sizeof(float),sizeof(double));
236 gmx_fio_setprecision(fio,bDouble);
237 fprintf(stderr,"Reading file %s, %s (%s precision)\n",
238 gmx_fio_getname(fio),buf,bDouble ? "double" : "single");
240 gmx_fio_do_int(fio,fver);
247 void set_inbox(int natom, rvec *x)
252 tmp[XX]=tmp[YY]=tmp[ZZ]=0.0;
255 if(x[i][XX]<tmp[XX]) tmp[XX]=x[i][XX];
256 if(x[i][YY]<tmp[YY]) tmp[YY]=x[i][YY];
257 if(x[i][ZZ]<tmp[ZZ]) tmp[ZZ]=x[i][ZZ];
264 int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist)
271 snew(tlist->index,at->nr);
272 for (i=0;i<at->nr;i++)
275 mol_id = get_mol_id(at->index[i],mtop->nmolblock,mtop->molblock,&type,&block);
278 if(tlist->index[j]==type)
283 tlist->index[nr]=type;
288 srenew(tlist->index,nr);
292 void check_types(t_block *ins_at,t_block *rest_at,gmx_mtop_t *mtop)
294 t_block *ins_mtype,*rest_mtype;
299 ins_mtype->nr = get_mtype_list(ins_at , mtop, ins_mtype );
300 rest_mtype->nr = get_mtype_list(rest_at, mtop, rest_mtype);
302 for(i=0;i<ins_mtype->nr;i++)
304 for(j=0;j<rest_mtype->nr;j++)
306 if(ins_mtype->index[i]==rest_mtype->index[j])
307 gmx_fatal(FARGS,"Moleculetype %s is found both in the group to insert and the rest of the system.\n"
308 "Because we need to exclude all interactions between the atoms in the group to\n"
309 "insert, the same moleculetype can not be used in both groups. Change the\n"
310 "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n"
311 "an appropriate *.itp file",*(mtop->moltype[rest_mtype->index[j]].name),
312 *(mtop->moltype[rest_mtype->index[j]].name));
316 sfree(ins_mtype->index);
317 sfree(rest_mtype->index);
322 int init_ins_at(t_block *ins_at,t_block *rest_at,t_state *state, pos_ins_t *pos_ins,gmx_groups_t *groups,int ins_grp_id, real xy_max)
325 real x,xmin,xmax,y,ymin,ymax,z,zmin,zmax;
327 snew(rest_at->index,state->natoms);
329 xmin=xmax=state->x[ins_at->index[0]][XX];
330 ymin=ymax=state->x[ins_at->index[0]][YY];
331 zmin=zmax=state->x[ins_at->index[0]][ZZ];
333 for(i=0;i<state->natoms;i++)
335 gid = groups->grpnr[egcFREEZE][i];
336 if(groups->grps[egcFREEZE].nm_ind[gid]==ins_grp_id)
354 srenew(rest_at->index,c);
358 pos_ins->xmin[XX]=xmin-((xmax-xmin)*xy_max-(xmax-xmin))/2;
359 pos_ins->xmin[YY]=ymin-((ymax-ymin)*xy_max-(ymax-ymin))/2;
361 pos_ins->xmax[XX]=xmax+((xmax-xmin)*xy_max-(xmax-xmin))/2;
362 pos_ins->xmax[YY]=ymax+((ymax-ymin)*xy_max-(ymax-ymin))/2;
364 pos_ins->xmin[XX]=xmin;
365 pos_ins->xmin[YY]=ymin;
367 pos_ins->xmax[XX]=xmax;
368 pos_ins->xmax[YY]=ymax;
371 /* 6.0 is estimated thickness of bilayer */
372 if( (zmax-zmin) < 6.0 )
374 pos_ins->xmin[ZZ]=zmin+(zmax-zmin)/2.0-3.0;
375 pos_ins->xmax[ZZ]=zmin+(zmax-zmin)/2.0+3.0;
377 pos_ins->xmin[ZZ]=zmin;
378 pos_ins->xmax[ZZ]=zmax;
384 real est_prot_area(pos_ins_t *pos_ins,rvec *r,t_block *ins_at, mem_t *mem_p)
386 real x,y,dx=0.15,dy=0.15,area=0.0;
390 for(x=pos_ins->xmin[XX];x<pos_ins->xmax[XX];x+=dx)
392 for(y=pos_ins->xmin[YY];y<pos_ins->xmax[YY];y+=dy)
399 if ( (r[at][XX]>=x) && (r[at][XX]<x+dx) &&
400 (r[at][YY]>=y) && (r[at][YY]<y+dy) &&
401 (r[at][ZZ]>mem_p->zmin+1.0) && (r[at][ZZ]<mem_p->zmax-1.0) )
404 } while ( (c<ins_at->nr) && (add<0.5) );
413 void init_lip(matrix box, gmx_mtop_t *mtop, lip_t *lip)
419 mem_area = box[XX][XX]*box[YY][YY]-box[XX][YY]*box[YY][XX];
420 for(i=0;i<mtop->nmolblock;i++)
422 if(mtop->molblock[i].type == lip->id)
424 lip->nr=mtop->molblock[i].nmol;
425 lip->natoms=mtop->molblock[i].natoms_mol;
428 lip->area=2.0*mem_area/(double)lip->nr;
430 for (i=0;i<lip->id;i++)
431 mol1+=mtop->molblock[i].nmol;
435 int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins)
437 int i,j,at,mol,nmol,nmolbox,count;
439 real z,zmin,zmax,mem_area;
445 mem_a=&(mem_p->mem_at);
446 snew(mol_id,mem_a->nr);
447 /* snew(index,mem_a->nr); */
448 zmin=pos_ins->xmax[ZZ];
449 zmax=pos_ins->xmin[ZZ];
450 for(i=0;i<mem_a->nr;i++)
453 if( (r[at][XX]>pos_ins->xmin[XX]) && (r[at][XX]<pos_ins->xmax[XX]) &&
454 (r[at][YY]>pos_ins->xmin[YY]) && (r[at][YY]<pos_ins->xmax[YY]) &&
455 (r[at][ZZ]>pos_ins->xmin[ZZ]) && (r[at][ZZ]<pos_ins->xmax[ZZ]) )
457 mol = get_mol_id(at,mtop->nmolblock,mtop->molblock,&type,&block);
474 /* index[count]=at;*/
481 mem_p->mol_id=mol_id;
482 /* srenew(index,count);*/
483 /* mem_p->mem_at.nr=count;*/
484 /* sfree(mem_p->mem_at.index);*/
485 /* mem_p->mem_at.index=index;*/
487 if((zmax-zmin)>(box[ZZ][ZZ]-0.5))
488 gmx_fatal(FARGS,"Something is wrong with your membrane. Max and min z values are %f and %f.\n"
489 "Maybe your membrane is not centered in the box, but located at the box edge in the z-direction,\n"
490 "so that one membrane is distributed over two periodic box images. Another possibility is that\n"
491 "your water layer is not thick enough.\n",zmax,zmin);
494 mem_p->zmed=(zmax-zmin)/2+zmin;
496 /*number of membrane molecules in protein box*/
497 nmolbox = count/mtop->molblock[block].natoms_mol;
498 /*mem_area = box[XX][XX]*box[YY][YY]-box[XX][YY]*box[YY][XX];
499 mem_p->lip_area = 2.0*mem_area/(double)mem_p->nmol;*/
500 mem_area = (pos_ins->xmax[XX]-pos_ins->xmin[XX])*(pos_ins->xmax[YY]-pos_ins->xmin[YY]);
501 mem_p->lip_area = 2.0*mem_area/(double)nmolbox;
503 return mem_p->mem_at.nr;
506 void init_resize(t_block *ins_at,rvec *r_ins,pos_ins_t *pos_ins,mem_t *mem_p,rvec *r, gmx_bool bALLOW_ASYMMETRY)
508 int i,j,at,c,outsidesum,gctr=0;
512 for (i=0;i<pos_ins->pieces;i++)
513 idxsum+=pos_ins->nidx[i];
514 if (idxsum!=ins_at->nr)
515 gmx_fatal(FARGS,"Piecewise sum of inserted atoms not same as size of group selected to insert.");
517 snew(pos_ins->geom_cent,pos_ins->pieces);
518 for (i=0;i<pos_ins->pieces;i++)
523 pos_ins->geom_cent[i][j]=0;
526 pos_ins->geom_cent[i][j]=0;
527 for (j=0;j<pos_ins->nidx[i];j++)
529 at=pos_ins->subindex[i][j];
530 copy_rvec(r[at],r_ins[gctr]);
531 if( (r_ins[gctr][ZZ]<mem_p->zmax) && (r_ins[gctr][ZZ]>mem_p->zmin) )
533 rvec_inc(pos_ins->geom_cent[i],r_ins[gctr]);
541 svmul(1/(double)c,pos_ins->geom_cent[i],pos_ins->geom_cent[i]);
542 if (!bALLOW_ASYMMETRY)
543 pos_ins->geom_cent[i][ZZ]=mem_p->zmed;
545 fprintf(stderr,"Embedding piece %d with center of geometry: %f %f %f\n",i,pos_ins->geom_cent[i][XX],pos_ins->geom_cent[i][YY],pos_ins->geom_cent[i][ZZ]);
547 fprintf(stderr,"\n");
550 void resize(t_block *ins_at, rvec *r_ins, rvec *r, pos_ins_t *pos_ins,rvec fac)
553 for (k=0;k<pos_ins->pieces;k++)
554 for(i=0;i<pos_ins->nidx[k];i++)
556 at=pos_ins->subindex[k][i];
558 r[at][j]=pos_ins->geom_cent[k][j]+fac[j]*(r_ins[c][j]-pos_ins->geom_cent[k][j]);
563 int gen_rm_list(rm_t *rm_p,t_block *ins_at,t_block *rest_at,t_pbc *pbc, gmx_mtop_t *mtop,
564 rvec *r, rvec *r_ins, mem_t *mem_p, pos_ins_t *pos_ins, real probe_rad, int low_up_rm, gmx_bool bALLOW_ASYMMETRY)
566 int i,j,k,l,at,at2,mol_id;
568 int nrm,nupper,nlower;
569 real r_min_rad,z_lip,min_norm;
575 r_min_rad=probe_rad*probe_rad;
576 snew(rm_p->mol,mtop->mols.nr);
577 snew(rm_p->block,mtop->mols.nr);
580 for(i=0;i<ins_at->nr;i++)
583 for(j=0;j<rest_at->nr;j++)
585 at2=rest_at->index[j];
586 pbc_dx(pbc,r[at],r[at2],dr);
588 if(norm2(dr)<r_min_rad)
590 mol_id = get_mol_id(at2,mtop->nmolblock,mtop->molblock,&type,&block);
593 if(rm_p->mol[l]==mol_id)
597 /*fprintf(stderr,"%d wordt toegevoegd\n",mol_id);*/
598 rm_p->mol[nrm]=mol_id;
599 rm_p->block[nrm]=block;
602 for(l=0;l<mem_p->nmol;l++)
604 if(mol_id==mem_p->mol_id[l])
606 for(k=mtop->mols.index[mol_id];k<mtop->mols.index[mol_id+1];k++)
608 z_lip/=mtop->molblock[block].natoms_mol;
609 if(z_lip<mem_p->zmed)
620 /*make sure equal number of lipids from upper and lower layer are removed */
621 if( (nupper!=nlower) && (!bALLOW_ASYMMETRY) )
623 snew(dist,mem_p->nmol);
624 snew(order,mem_p->nmol);
625 for(i=0;i<mem_p->nmol;i++)
627 at = mtop->mols.index[mem_p->mol_id[i]];
628 pbc_dx(pbc,r[at],pos_ins->geom_cent[0],dr);
629 if (pos_ins->pieces>1)
633 for (k=1;k<pos_ins->pieces;k++)
635 pbc_dx(pbc,r[at],pos_ins->geom_cent[k],dr_tmp);
636 if (norm2(dr_tmp) < min_norm)
638 min_norm=norm2(dr_tmp);
639 copy_rvec(dr_tmp,dr);
643 dist[i]=dr[XX]*dr[XX]+dr[YY]*dr[YY];
645 while (j>=0 && dist[i]<dist[order[j]])
654 while(nupper!=nlower)
656 mol_id=mem_p->mol_id[order[i]];
657 block=get_block(mol_id,mtop->nmolblock,mtop->molblock);
661 if(rm_p->mol[l]==mol_id)
666 for(k=mtop->mols.index[mol_id];k<mtop->mols.index[mol_id+1];k++)
668 z_lip/=mtop->molblock[block].natoms_mol;
669 if(nupper>nlower && z_lip<mem_p->zmed)
671 rm_p->mol[nrm]=mol_id;
672 rm_p->block[nrm]=block;
676 else if (nupper<nlower && z_lip>mem_p->zmed)
678 rm_p->mol[nrm]=mol_id;
679 rm_p->block[nrm]=block;
687 gmx_fatal(FARGS,"Trying to remove more lipid molecules than there are in the membrane");
694 srenew(rm_p->mol,nrm);
695 srenew(rm_p->block,nrm);
697 return nupper+nlower;
700 void rm_group(t_inputrec *ir, gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state *state, t_block *ins_at, pos_ins_t *pos_ins)
702 int i,j,k,n,rm,mol_id,at,block;
704 atom_id *list,*new_mols;
705 unsigned char *new_egrp[egcNR];
708 snew(list,state->natoms);
710 for(i=0;i<rm_p->nr;i++)
713 at=mtop->mols.index[mol_id];
714 block =rm_p->block[i];
715 mtop->molblock[block].nmol--;
716 for(j=0;j<mtop->molblock[block].natoms_mol;j++)
722 mtop->mols.index[mol_id]=-1;
725 mtop->mols.nr-=rm_p->nr;
726 mtop->mols.nalloc_index-=rm_p->nr;
727 snew(new_mols,mtop->mols.nr);
728 for(i=0;i<mtop->mols.nr+rm_p->nr;i++)
731 if(mtop->mols.index[i]!=-1)
733 new_mols[j]=mtop->mols.index[i];
737 sfree(mtop->mols.index);
738 mtop->mols.index=new_mols;
743 state->nalloc=state->natoms;
744 snew(x_tmp,state->nalloc);
745 snew(v_tmp,state->nalloc);
749 if(groups->grpnr[i]!=NULL)
751 groups->ngrpnr[i]=state->natoms;
752 snew(new_egrp[i],state->natoms);
757 for (i=0;i<state->natoms+n;i++)
773 if(groups->grpnr[j]!=NULL)
775 new_egrp[j][i-rm]=groups->grpnr[j][i];
778 copy_rvec(state->x[i],x_tmp[i-rm]);
779 copy_rvec(state->v[i],v_tmp[i-rm]);
780 for(j=0;j<ins_at->nr;j++)
782 if (i==ins_at->index[j])
783 ins_at->index[j]=i-rm;
785 for(j=0;j<pos_ins->pieces;j++)
787 for(k=0;k<pos_ins->nidx[j];k++)
789 if (i==pos_ins->subindex[j][k])
790 pos_ins->subindex[j][k]=i-rm;
802 if(groups->grpnr[i]!=NULL)
804 sfree(groups->grpnr[i]);
805 groups->grpnr[i]=new_egrp[i];
810 int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop)
813 int type,natom,nmol,at,atom1=0,rm_at=0;
815 /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/
816 /*this routine does not live as dangerously as it seems. There is namely a check in mdrunner_membed to make
817 *sure that g_membed exits with a warning when there are molecules of the same type not in the
818 *ins_at index group. MGWolf 050710 */
821 snew(bRM,mtop->nmoltype);
822 for (i=0;i<mtop->nmoltype;i++)
827 for (i=0;i<mtop->nmolblock;i++)
829 /*loop over molecule blocks*/
830 type =mtop->molblock[i].type;
831 natom =mtop->molblock[i].natoms_mol;
832 nmol =mtop->molblock[i].nmol;
834 for(j=0;j<natom*nmol && bRM[type]==TRUE;j++)
836 /*loop over atoms in the block*/
837 at=j+atom1; /*atom index = block index + offset*/
840 for (m=0;(m<ins_at->nr) && (bINS==FALSE);m++)
842 /*loop over atoms in insertion index group to determine if we're inserting one*/
843 if(at==ins_at->index[m])
850 atom1+=natom*nmol; /*update offset*/
853 rm_at+=natom*nmol; /*increment bonded removal counter by # atoms in block*/
857 for(i=0;i<mtop->nmoltype;i++)
863 mtop->moltype[i].ilist[j].nr=0;
865 for(j=F_POSRES;j<=F_VSITEN;j++)
867 mtop->moltype[i].ilist[j].nr=0;
876 void top_update(const char *topfile, char *ins, rm_t *rm_p, gmx_mtop_t *mtop)
878 #define TEMP_FILENM "temp.top"
881 char buf[STRLEN],buf2[STRLEN],*temp;
882 int i,*nmol_rm,nmol,line;
884 fpin = ffopen(topfile,"r");
885 fpout = ffopen(TEMP_FILENM,"w");
887 snew(nmol_rm,mtop->nmoltype);
888 for(i=0;i<rm_p->nr;i++)
889 nmol_rm[rm_p->block[i]]++;
892 while(fgets(buf,STRLEN,fpin))
898 if ((temp=strchr(buf2,'\n')) != NULL)
905 if ((temp=strchr(buf2,'\n')) != NULL)
908 if (buf2[strlen(buf2)-1]==']')
910 buf2[strlen(buf2)-1]='\0';
913 if (gmx_strcasecmp(buf2,"molecules")==0)
916 fprintf(fpout,"%s",buf);
917 } else if (bMolecules==1)
919 for(i=0;i<mtop->nmolblock;i++)
921 nmol=mtop->molblock[i].nmol;
922 sprintf(buf,"%-15s %5d\n",*(mtop->moltype[mtop->molblock[i].type].name),nmol);
923 fprintf(fpout,"%s",buf);
926 } else if (bMolecules==2)
931 fprintf(fpout,"%s",buf);
935 fprintf(fpout,"%s",buf);
940 /* use ffopen to generate backup of topinout */
941 fpout=ffopen(topfile,"w");
943 rename(TEMP_FILENM,topfile);
947 void md_print_warning(const t_commrec *cr,FILE *fplog,const char *buf)
951 fprintf(stderr,"\n%s\n",buf);
955 fprintf(fplog,"\n%s\n",buf);
959 /* simulation conditions to transmit. Keep in mind that they are
960 transmitted to other nodes through an MPI_Reduce after
961 casting them to a real (so the signals can be sent together with other
962 data). This means that the only meaningful values are positive,
964 enum { eglsNABNSB, eglsCHKPT, eglsSTOPCOND, eglsRESETCOUNTERS, eglsNR };
965 /* Is the signal in one simulation independent of other simulations? */
966 gmx_bool gs_simlocal[eglsNR] = { TRUE, FALSE, FALSE, TRUE };
969 int nstms; /* The frequency for intersimulation communication */
970 int sig[eglsNR]; /* The signal set by one process in do_md */
971 int set[eglsNR]; /* The communicated signal, equal for all processes */
975 static int multisim_min(const gmx_multisim_t *ms,int nmin,int n)
978 gmx_bool bPos,bEqual;
983 gmx_sumi_sim(ms->nsim,buf,ms);
986 for(s=0; s<ms->nsim; s++)
988 bPos = bPos && (buf[s] > 0);
989 bEqual = bEqual && (buf[s] == buf[0]);
995 nmin = min(nmin,buf[0]);
999 /* Find the least common multiple */
1000 for(d=2; d<nmin; d++)
1003 while (s < ms->nsim && d % buf[s] == 0)
1009 /* We found the LCM and it is less than nmin */
1021 static int multisim_nstsimsync(const t_commrec *cr,
1022 const t_inputrec *ir,int repl_ex_nst)
1029 nmin = multisim_min(cr->ms,nmin,ir->nstlist);
1030 nmin = multisim_min(cr->ms,nmin,ir->nstcalcenergy);
1031 nmin = multisim_min(cr->ms,nmin,repl_ex_nst);
1032 if (nmin == INT_MAX)
1034 gmx_fatal(FARGS,"Can not find an appropriate interval for inter-simulation communication, since nstlist, nstcalcenergy and -replex are all <= 0");
1036 /* Avoid inter-simulation communication at every (second) step */
1043 gmx_bcast(sizeof(int),&nmin,cr);
1048 static void init_global_signals(globsig_t *gs,const t_commrec *cr,
1049 const t_inputrec *ir,int repl_ex_nst)
1055 gs->nstms = multisim_nstsimsync(cr,ir,repl_ex_nst);
1058 fprintf(debug,"Syncing simulations for checkpointing and termination every %d steps\n",gs->nstms);
1066 for(i=0; i<eglsNR; i++)
1073 static void copy_coupling_state(t_state *statea,t_state *stateb,
1074 gmx_ekindata_t *ekinda,gmx_ekindata_t *ekindb, t_grpopts* opts)
1077 /* MRS note -- might be able to get rid of some of the arguments. Look over it when it's all debugged */
1081 /* Make sure we have enough space for x and v */
1082 if (statea->nalloc > stateb->nalloc)
1084 stateb->nalloc = statea->nalloc;
1085 srenew(stateb->x,stateb->nalloc);
1086 srenew(stateb->v,stateb->nalloc);
1089 stateb->natoms = statea->natoms;
1090 stateb->ngtc = statea->ngtc;
1091 stateb->nnhpres = statea->nnhpres;
1092 stateb->veta = statea->veta;
1095 copy_mat(ekinda->ekin,ekindb->ekin);
1096 for (i=0; i<stateb->ngtc; i++)
1098 ekindb->tcstat[i].T = ekinda->tcstat[i].T;
1099 ekindb->tcstat[i].Th = ekinda->tcstat[i].Th;
1100 copy_mat(ekinda->tcstat[i].ekinh,ekindb->tcstat[i].ekinh);
1101 copy_mat(ekinda->tcstat[i].ekinf,ekindb->tcstat[i].ekinf);
1102 ekindb->tcstat[i].ekinscalef_nhc = ekinda->tcstat[i].ekinscalef_nhc;
1103 ekindb->tcstat[i].ekinscaleh_nhc = ekinda->tcstat[i].ekinscaleh_nhc;
1104 ekindb->tcstat[i].vscale_nhc = ekinda->tcstat[i].vscale_nhc;
1107 copy_rvecn(statea->x,stateb->x,0,stateb->natoms);
1108 copy_rvecn(statea->v,stateb->v,0,stateb->natoms);
1109 copy_mat(statea->box,stateb->box);
1110 copy_mat(statea->box_rel,stateb->box_rel);
1111 copy_mat(statea->boxv,stateb->boxv);
1113 for (i = 0; i<stateb->ngtc; i++)
1115 nc = i*opts->nhchainlength;
1116 for (j=0; j<opts->nhchainlength; j++)
1118 stateb->nosehoover_xi[nc+j] = statea->nosehoover_xi[nc+j];
1119 stateb->nosehoover_vxi[nc+j] = statea->nosehoover_vxi[nc+j];
1122 if (stateb->nhpres_xi != NULL)
1124 for (i = 0; i<stateb->nnhpres; i++)
1126 nc = i*opts->nhchainlength;
1127 for (j=0; j<opts->nhchainlength; j++)
1129 stateb->nhpres_xi[nc+j] = statea->nhpres_xi[nc+j];
1130 stateb->nhpres_vxi[nc+j] = statea->nhpres_vxi[nc+j];
1136 static void compute_globals(FILE *fplog, gmx_global_stat_t gstat, t_commrec *cr, t_inputrec *ir,
1137 t_forcerec *fr, gmx_ekindata_t *ekind,
1138 t_state *state, t_state *state_global, t_mdatoms *mdatoms,
1139 t_nrnb *nrnb, t_vcm *vcm, gmx_wallcycle_t wcycle,
1140 gmx_enerdata_t *enerd,tensor force_vir, tensor shake_vir, tensor total_vir,
1141 tensor pres, rvec mu_tot, gmx_constr_t constr,
1142 globsig_t *gs,gmx_bool bInterSimGS,
1143 matrix box, gmx_mtop_t *top_global, real *pcurr,
1144 int natoms, gmx_bool *bSumEkinhOld, int flags)
1147 real gs_buf[eglsNR];
1148 tensor corr_vir,corr_pres;
1149 gmx_bool bEner,bPres,bTemp;
1150 gmx_bool bRerunMD, bStopCM, bGStat, bIterate,
1151 bFirstIterate,bReadEkin,bEkinAveVel,bScaleEkin, bConstrain;
1152 real prescorr,enercorr,dvdlcorr;
1154 /* translate CGLO flags to gmx_booleans */
1155 bRerunMD = flags & CGLO_RERUNMD;
1156 bStopCM = flags & CGLO_STOPCM;
1157 bGStat = flags & CGLO_GSTAT;
1158 bReadEkin = (flags & CGLO_READEKIN);
1159 bScaleEkin = (flags & CGLO_SCALEEKIN);
1160 bEner = flags & CGLO_ENERGY;
1161 bTemp = flags & CGLO_TEMPERATURE;
1162 bPres = (flags & CGLO_PRESSURE);
1163 bConstrain = (flags & CGLO_CONSTRAINT);
1164 bIterate = (flags & CGLO_ITERATE);
1165 bFirstIterate = (flags & CGLO_FIRSTITERATE);
1167 /* we calculate a full state kinetic energy either with full-step velocity verlet
1168 or half step where we need the pressure */
1169 bEkinAveVel = (ir->eI==eiVV || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir) && bPres) || bReadEkin);
1171 /* in initalization, it sums the shake virial in vv, and to
1172 sums ekinh_old in leapfrog (or if we are calculating ekinh_old for other reasons */
1174 /* ########## Kinetic energy ############## */
1178 /* Non-equilibrium MD: this is parallellized, but only does communication
1179 * when there really is NEMD.
1182 if (PAR(cr) && (ekind->bNEMD))
1184 accumulate_u(cr,&(ir->opts),ekind);
1189 restore_ekinstate_from_state(cr,ekind,&state_global->ekinstate);
1194 calc_ke_part(state,&(ir->opts),mdatoms,ekind,nrnb,bEkinAveVel,bIterate);
1199 /* Calculate center of mass velocity if necessary, also parallellized */
1200 if (bStopCM && !bRerunMD && bEner)
1202 calc_vcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms,
1203 state->x,state->v,vcm);
1207 if (bTemp || bPres || bEner || bConstrain)
1211 /* We will not sum ekinh_old,
1212 * so signal that we still have to do it.
1214 *bSumEkinhOld = TRUE;
1221 for(i=0; i<eglsNR; i++)
1223 gs_buf[i] = gs->sig[i];
1228 wallcycle_start(wcycle,ewcMoveE);
1229 GMX_MPE_LOG(ev_global_stat_start);
1230 global_stat(fplog,gstat,cr,enerd,force_vir,shake_vir,mu_tot,
1231 ir,ekind,constr,vcm,
1232 gs != NULL ? eglsNR : 0,gs_buf,
1234 *bSumEkinhOld,flags);
1235 GMX_MPE_LOG(ev_global_stat_finish);
1236 wallcycle_stop(wcycle,ewcMoveE);
1240 if (MULTISIM(cr) && bInterSimGS)
1244 /* Communicate the signals between the simulations */
1245 gmx_sum_sim(eglsNR,gs_buf,cr->ms);
1247 /* Communicate the signals form the master to the others */
1248 gmx_bcast(eglsNR*sizeof(gs_buf[0]),gs_buf,cr);
1250 for(i=0; i<eglsNR; i++)
1252 if (bInterSimGS || gs_simlocal[i])
1254 /* Set the communicated signal only when it is non-zero,
1255 * since signals might not be processed at each MD step.
1257 gsi = (gs_buf[i] >= 0 ?
1258 (int)(gs_buf[i] + 0.5) :
1259 (int)(gs_buf[i] - 0.5));
1264 /* Turn off the local signal */
1269 *bSumEkinhOld = FALSE;
1273 if (!ekind->bNEMD && debug && bTemp && (vcm->nr > 0))
1276 mdatoms->start,mdatoms->start+mdatoms->homenr,
1277 state->v,vcm->group_p[0],
1278 mdatoms->massT,mdatoms->tmass,ekind->ekin);
1282 /* Do center of mass motion removal */
1283 if (bStopCM && !bRerunMD) /* is this correct? Does it get called too often with this logic? */
1285 check_cm_grp(fplog,vcm,ir,1);
1286 do_stopcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms->cVCM,
1287 state->x,state->v,vcm);
1288 inc_nrnb(nrnb,eNR_STOPCM,mdatoms->homenr);
1294 /* Sum the kinetic energies of the groups & calc temp */
1295 /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with IR_NPT_TROTTER */
1296 /* three maincase: VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
1297 Leap with AveVel is also an option for the future but not supported now.
1298 bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
1299 If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
1300 bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
1301 If FALSE, we go ahead and erase over it.
1303 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,&(enerd->term[F_DKDL]),
1304 bEkinAveVel,bIterate,bScaleEkin);
1306 enerd->term[F_EKIN] = trace(ekind->ekin);
1309 /* ########## Long range energy information ###### */
1311 if (bEner || bPres || bConstrain)
1313 calc_dispcorr(fplog,ir,fr,0,top_global->natoms,box,state->lambda,
1314 corr_pres,corr_vir,&prescorr,&enercorr,&dvdlcorr);
1317 if (bEner && bFirstIterate)
1319 enerd->term[F_DISPCORR] = enercorr;
1320 enerd->term[F_EPOT] += enercorr;
1321 enerd->term[F_DVDL] += dvdlcorr;
1322 if (fr->efep != efepNO) {
1323 enerd->dvdl_lin += dvdlcorr;
1327 /* ########## Now pressure ############## */
1328 if (bPres || bConstrain)
1331 m_add(force_vir,shake_vir,total_vir);
1333 /* Calculate pressure and apply LR correction if PPPM is used.
1334 * Use the box from last timestep since we already called update().
1337 enerd->term[F_PRES] = calc_pres(fr->ePBC,ir->nwall,box,ekind->ekin,total_vir,pres,
1338 (fr->eeltype==eelPPPM)?enerd->term[F_COUL_RECIP]:0.0);
1340 /* Calculate long range corrections to pressure and energy */
1341 /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
1342 and computes enerd->term[F_DISPCORR]. Also modifies the
1343 total_vir and pres tesors */
1345 m_add(total_vir,corr_vir,total_vir);
1346 m_add(pres,corr_pres,pres);
1347 enerd->term[F_PDISPCORR] = prescorr;
1348 enerd->term[F_PRES] += prescorr;
1349 *pcurr = enerd->term[F_PRES];
1350 /* calculate temperature using virial */
1351 enerd->term[F_VTEMP] = calc_temp(trace(total_vir),ir->opts.nrdf[0]);
1357 /* Definitions for convergence of iterated constraints */
1359 /* iterate constraints up to 50 times */
1360 #define MAXITERCONST 50
1365 real f,fprev,x,xprev;
1368 real allrelerr[MAXITERCONST+2];
1369 int num_close; /* number of "close" violations, caused by limited precision. */
1373 #define CONVERGEITER 0.000000001
1374 #define CLOSE_ENOUGH 0.000001000
1376 #define CONVERGEITER 0.0001
1377 #define CLOSE_ENOUGH 0.0050
1380 /* we want to keep track of the close calls. If there are too many, there might be some other issues.
1381 so we make sure that it's either less than some predetermined number, or if more than that number,
1382 only some small fraction of the total. */
1383 #define MAX_NUMBER_CLOSE 50
1384 #define FRACTION_CLOSE 0.001
1386 /* maximum length of cyclic traps to check, emerging from limited numerical precision */
1389 static void gmx_iterate_init(gmx_iterate_t *iterate,gmx_bool bIterate)
1393 iterate->iter_i = 0;
1394 iterate->bIterate = bIterate;
1395 iterate->num_close = 0;
1396 for (i=0;i<MAXITERCONST+2;i++)
1398 iterate->allrelerr[i] = 0;
1402 static gmx_bool done_iterating(const t_commrec *cr,FILE *fplog, int nsteps, gmx_iterate_t *iterate, gmx_bool bFirstIterate, real fom, real *newf)
1404 /* monitor convergence, and use a secant search to propose new
1407 The secant method computes x_{i+1} = x_{i} - f(x_{i}) * ---------------------
1408 f(x_{i}) - f(x_{i-1})
1410 The function we are trying to zero is fom-x, where fom is the
1411 "figure of merit" which is the pressure (or the veta value) we
1412 would get by putting in an old value of the pressure or veta into
1413 the incrementor function for the step or half step. I have
1414 verified that this gives the same answer as self consistent
1415 iteration, usually in many fewer steps, especially for small tau_p.
1417 We could possibly eliminate an iteration with proper use
1418 of the value from the previous step, but that would take a bit
1419 more bookkeeping, especially for veta, since tests indicate the
1420 function of veta on the last step is not sufficiently close to
1421 guarantee convergence this step. This is
1422 good enough for now. On my tests, I could use tau_p down to
1423 0.02, which is smaller that would ever be necessary in
1424 practice. Generally, 3-5 iterations will be sufficient */
1434 iterate->f = fom-iterate->x;
1441 iterate->f = fom-iterate->x; /* we want to zero this difference */
1442 if ((iterate->iter_i > 1) && (iterate->iter_i < MAXITERCONST))
1444 if (iterate->f==iterate->fprev)
1450 *newf = iterate->x - (iterate->x-iterate->xprev)*(iterate->f)/(iterate->f-iterate->fprev);
1455 /* just use self-consistent iteration the first step to initialize, or
1456 if it's not converging (which happens occasionally -- need to investigate why) */
1460 /* Consider a slight shortcut allowing us to exit one sooner -- we check the
1461 difference between the closest of x and xprev to the new
1462 value. To be 100% certain, we should check the difference between
1463 the last result, and the previous result, or
1465 relerr = (fabs((x-xprev)/fom));
1467 but this is pretty much never necessary under typical conditions.
1468 Checking numerically, it seems to lead to almost exactly the same
1469 trajectories, but there are small differences out a few decimal
1470 places in the pressure, and eventually in the v_eta, but it could
1473 if (fabs(*newf-x) < fabs(*newf - xprev)) { xmin = x;} else { xmin = xprev;}
1474 relerr = (fabs((*newf-xmin) / *newf));
1477 err = fabs((iterate->f-iterate->fprev));
1478 relerr = fabs(err/fom);
1480 iterate->allrelerr[iterate->iter_i] = relerr;
1482 if (iterate->iter_i > 0)
1486 fprintf(debug,"Iterating NPT constraints: %6i %20.12f%14.6g%20.12f\n",
1487 iterate->iter_i,fom,relerr,*newf);
1490 if ((relerr < CONVERGEITER) || (err < CONVERGEITER) || (fom==0) || ((iterate->x == iterate->xprev) && iterate->iter_i > 1))
1492 iterate->bIterate = FALSE;
1495 fprintf(debug,"Iterating NPT constraints: CONVERGED\n");
1499 if (iterate->iter_i > MAXITERCONST)
1501 if (relerr < CLOSE_ENOUGH)
1504 for (i=1;i<CYCLEMAX;i++) {
1505 if ((iterate->allrelerr[iterate->iter_i-(1+i)] == iterate->allrelerr[iterate->iter_i-1]) &&
1506 (iterate->allrelerr[iterate->iter_i-(1+i)] == iterate->allrelerr[iterate->iter_i-(1+2*i)])) {
1510 fprintf(debug,"Exiting from an NPT iterating cycle of length %d\n",i);
1517 /* step 1: trapped in a numerical attractor */
1518 /* we are trapped in a numerical attractor, and can't converge any more, and are close to the final result.
1519 Better to give up convergence here than have the simulation die.
1521 iterate->num_close++;
1526 /* Step #2: test if we are reasonably close for other reasons, then monitor the number. If not, die */
1528 /* how many close calls have we had? If less than a few, we're OK */
1529 if (iterate->num_close < MAX_NUMBER_CLOSE)
1531 sprintf(buf,"Slight numerical convergence deviation with NPT at step %d, relative error only %10.5g, likely not a problem, continuing\n",nsteps,relerr);
1532 md_print_warning(cr,fplog,buf);
1533 iterate->num_close++;
1535 /* if more than a few, check the total fraction. If too high, die. */
1536 } else if (iterate->num_close/(double)nsteps > FRACTION_CLOSE) {
1537 gmx_fatal(FARGS,"Could not converge NPT constraints, too many exceptions (%d%%\n",iterate->num_close/(double)nsteps);
1543 gmx_fatal(FARGS,"Could not converge NPT constraints\n");
1548 iterate->xprev = iterate->x;
1550 iterate->fprev = iterate->f;
1556 static void check_nst_param(FILE *fplog,t_commrec *cr,
1557 const char *desc_nst,int nst,
1558 const char *desc_p,int *p)
1562 if (*p > 0 && *p % nst != 0)
1564 /* Round up to the next multiple of nst */
1565 *p = ((*p)/nst + 1)*nst;
1566 sprintf(buf,"NOTE: %s changes %s to %d\n",desc_nst,desc_p,*p);
1567 md_print_warning(cr,fplog,buf);
1571 static void reset_all_counters(FILE *fplog,t_commrec *cr,
1572 gmx_large_int_t step,
1573 gmx_large_int_t *step_rel,t_inputrec *ir,
1574 gmx_wallcycle_t wcycle,t_nrnb *nrnb,
1575 gmx_runtime_t *runtime)
1577 char buf[STRLEN],sbuf[STEPSTRSIZE];
1579 /* Reset all the counters related to performance over the run */
1580 sprintf(buf,"Step %s: resetting all time and cycle counters\n",
1581 gmx_step_str(step,sbuf));
1582 md_print_warning(cr,fplog,buf);
1584 wallcycle_stop(wcycle,ewcRUN);
1585 wallcycle_reset_all(wcycle);
1586 if (DOMAINDECOMP(cr))
1588 reset_dd_statistics_counters(cr->dd);
1591 ir->init_step += *step_rel;
1592 ir->nsteps -= *step_rel;
1594 wallcycle_start(wcycle,ewcRUN);
1595 runtime_start(runtime);
1596 print_date_and_time(fplog,cr->nodeid,"Restarted time",runtime);
1599 static int check_nstglobalcomm(FILE *fplog,t_commrec *cr,
1600 int nstglobalcomm,t_inputrec *ir)
1604 if (!EI_DYNAMICS(ir->eI))
1609 if (nstglobalcomm == -1)
1611 if (ir->nstcalcenergy == 0 && ir->nstlist == 0)
1614 if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
1616 nstglobalcomm = ir->nstenergy;
1621 /* We assume that if nstcalcenergy > nstlist,
1622 * nstcalcenergy is a multiple of nstlist.
1624 if (ir->nstcalcenergy == 0 ||
1625 (ir->nstlist > 0 && ir->nstlist < ir->nstcalcenergy))
1627 nstglobalcomm = ir->nstlist;
1631 nstglobalcomm = ir->nstcalcenergy;
1637 if (ir->nstlist > 0 &&
1638 nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
1640 nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
1641 sprintf(buf,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm);
1642 md_print_warning(cr,fplog,buf);
1644 if (nstglobalcomm > ir->nstcalcenergy)
1646 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
1647 "nstcalcenergy",&ir->nstcalcenergy);
1650 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
1651 "nstenergy",&ir->nstenergy);
1653 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
1654 "nstlog",&ir->nstlog);
1657 if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
1659 sprintf(buf,"WARNING: Changing nstcomm from %d to %d\n",
1660 ir->nstcomm,nstglobalcomm);
1661 md_print_warning(cr,fplog,buf);
1662 ir->nstcomm = nstglobalcomm;
1665 return nstglobalcomm;
1668 void check_ir_old_tpx_versions(t_commrec *cr,FILE *fplog,
1669 t_inputrec *ir,gmx_mtop_t *mtop)
1671 /* Check required for old tpx files */
1672 if (IR_TWINRANGE(*ir) && ir->nstlist > 1 &&
1673 ir->nstcalcenergy % ir->nstlist != 0)
1675 md_print_warning(cr,fplog,"Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies");
1677 if (gmx_mtop_ftype_count(mtop,F_CONSTR) +
1678 gmx_mtop_ftype_count(mtop,F_CONSTRNC) > 0 &&
1679 ir->eConstrAlg == econtSHAKE)
1681 md_print_warning(cr,fplog,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect");
1682 if (ir->epc != epcNO)
1684 gmx_fatal(FARGS,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
1687 check_nst_param(fplog,cr,"nstlist",ir->nstlist,
1688 "nstcalcenergy",&ir->nstcalcenergy);
1689 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
1690 "nstenergy",&ir->nstenergy);
1691 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
1692 "nstlog",&ir->nstlog);
1693 if (ir->efep != efepNO)
1695 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
1696 "nstdhdl",&ir->nstdhdl);
1702 gmx_bool bGStatEveryStep;
1703 gmx_large_int_t step_ns;
1704 gmx_large_int_t step_nscheck;
1705 gmx_large_int_t nns;
1715 static void reset_nlistheuristics(gmx_nlheur_t *nlh,gmx_large_int_t step)
1719 nlh->step_nscheck = step;
1722 static void init_nlistheuristics(gmx_nlheur_t *nlh,
1723 gmx_bool bGStatEveryStep,gmx_large_int_t step)
1725 nlh->bGStatEveryStep = bGStatEveryStep;
1732 reset_nlistheuristics(nlh,step);
1735 static void update_nliststatistics(gmx_nlheur_t *nlh,gmx_large_int_t step)
1737 gmx_large_int_t nl_lt;
1738 char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
1740 /* Determine the neighbor list life time */
1741 nl_lt = step - nlh->step_ns;
1744 fprintf(debug,"%d atoms beyond ns buffer, updating neighbor list after %s steps\n",nlh->nabnsb,gmx_step_str(nl_lt,sbuf));
1748 nlh->s2 += nl_lt*nl_lt;
1749 nlh->ab += nlh->nabnsb;
1750 if (nlh->lt_runav == 0)
1752 nlh->lt_runav = nl_lt;
1753 /* Initialize the fluctuation average
1754 * such that at startup we check after 0 steps.
1756 nlh->lt_runav2 = sqr(nl_lt/2.0);
1758 /* Running average with 0.9 gives an exp. history of 9.5 */
1759 nlh->lt_runav2 = 0.9*nlh->lt_runav2 + 0.1*sqr(nlh->lt_runav - nl_lt);
1760 nlh->lt_runav = 0.9*nlh->lt_runav + 0.1*nl_lt;
1761 if (nlh->bGStatEveryStep)
1763 /* Always check the nlist validity */
1764 nlh->step_nscheck = step;
1768 /* We check after: <life time> - 2*sigma
1769 * The factor 2 is quite conservative,
1770 * but we assume that with nstlist=-1 the user
1771 * prefers exact integration over performance.
1773 nlh->step_nscheck = step
1774 + (int)(nlh->lt_runav - 2.0*sqrt(nlh->lt_runav2)) - 1;
1778 fprintf(debug,"nlist life time %s run av. %4.1f sig %3.1f check %s check with -gcom %d\n",
1779 gmx_step_str(nl_lt,sbuf),nlh->lt_runav,sqrt(nlh->lt_runav2),
1780 gmx_step_str(nlh->step_nscheck-step+1,sbuf2),
1781 (int)(nlh->lt_runav - 2.0*sqrt(nlh->lt_runav2)));
1785 static void set_nlistheuristics(gmx_nlheur_t *nlh,gmx_bool bReset,gmx_large_int_t step)
1791 reset_nlistheuristics(nlh,step);
1795 update_nliststatistics(nlh,step);
1798 nlh->step_ns = step;
1799 /* Initialize the cumulative coordinate scaling matrix */
1800 clear_mat(nlh->scale_tot);
1801 for(d=0; d<DIM; d++)
1803 nlh->scale_tot[d][d] = 1.0;
1807 double do_md_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
1808 const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
1810 gmx_vsite_t *vsite,gmx_constr_t constr,
1811 int stepout,t_inputrec *ir,
1812 gmx_mtop_t *top_global,
1814 t_state *state_global,
1816 t_nrnb *nrnb,gmx_wallcycle_t wcycle,
1817 gmx_edsam_t ed,t_forcerec *fr,
1818 int repl_ex_nst,int repl_ex_seed,
1819 real cpt_period,real max_hours,
1820 const char *deviceOptions,
1821 unsigned long Flags,
1822 gmx_runtime_t *runtime,
1823 rvec fac, rvec *r_ins, pos_ins_t *pos_ins, t_block *ins_at,
1824 real xy_step, real z_step, int it_xy, int it_z)
1827 gmx_large_int_t step,step_rel;
1830 gmx_bool bGStatEveryStep,bGStat,bNstEner,bCalcEnerPres;
1831 gmx_bool bNS,bNStList,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE,
1832 bFirstStep,bStateFromTPX,bInitStep,bLastStep,
1833 bBornRadii,bStartingFromCpt;
1834 gmx_bool bDoDHDL=FALSE;
1835 gmx_bool do_ene,do_log,do_verbose,bRerunWarnNoV=TRUE,
1836 bForceUpdate=FALSE,bCPT;
1838 gmx_bool bMasterState;
1839 int force_flags,cglo_flags;
1840 tensor force_vir,shake_vir,total_vir,tmp_vir,pres;
1842 t_trxstatus *status;
1845 t_state *bufstate=NULL;
1846 matrix *scale_tot,pcoupl_mu,M,ebox;
1848 t_trxframe rerun_fr;
1849 /* gmx_repl_ex_t repl_ex=NULL;*/
1852 gmx_localtop_t *top;
1853 t_mdebin *mdebin=NULL;
1854 t_state *state=NULL;
1855 rvec *f_global=NULL;
1858 gmx_enerdata_t *enerd;
1860 gmx_global_stat_t gstat;
1861 gmx_update_t upd=NULL;
1862 t_graph *graph=NULL;
1866 gmx_groups_t *groups;
1867 gmx_ekindata_t *ekind, *ekind_save;
1868 gmx_shellfc_t shellfc;
1869 int count,nconverged=0;
1872 gmx_bool bIonize=FALSE;
1873 gmx_bool bTCR=FALSE,bConverged=TRUE,bOK,bSumEkinhOld,bExchanged;
1875 gmx_bool bResetCountersHalfMaxH=FALSE;
1876 gmx_bool bVV,bIterations,bIterate,bFirstIterate,bTemp,bPres,bTrotter;
1879 rvec *xcopy=NULL,*vcopy=NULL,*cbuf=NULL;
1880 matrix boxcopy={{0}},lastbox;
1881 real veta_save,pcurr,scalevir,tracevir;
1884 real last_conserved = 0;
1888 char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
1889 int handled_stop_condition=gmx_stop_cond_none; /* compare to get_stop_condition*/
1890 gmx_iterate_t iterate;
1892 /* Temporary addition for FAHCORE checkpointing */
1896 /* Check for special mdrun options */
1897 bRerunMD = (Flags & MD_RERUN);
1898 bIonize = (Flags & MD_IONIZE);
1899 bFFscan = (Flags & MD_FFSCAN);
1900 bAppend = (Flags & MD_APPENDFILES);
1901 bGStatEveryStep = FALSE;
1902 if (Flags & MD_RESETCOUNTERSHALFWAY)
1906 /* Signal to reset the counters half the simulation steps. */
1907 wcycle_set_reset_counters(wcycle,ir->nsteps/2);
1909 /* Signal to reset the counters halfway the simulation time. */
1910 bResetCountersHalfMaxH = (max_hours > 0);
1913 /* md-vv uses averaged full step velocities for T-control
1914 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
1915 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
1916 bVV = EI_VV(ir->eI);
1917 if (bVV) /* to store the initial velocities while computing virial */
1919 snew(cbuf,top_global->natoms);
1921 /* all the iteratative cases - only if there are constraints */
1922 bIterations = ((IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
1923 bTrotter = (bVV && (IR_NPT_TROTTER(ir) || (IR_NVT_TROTTER(ir))));
1927 /* Since we don't know if the frames read are related in any way,
1928 * rebuild the neighborlist at every step.
1931 ir->nstcalcenergy = 1;
1935 check_ir_old_tpx_versions(cr,fplog,ir,top_global);
1937 nstglobalcomm = check_nstglobalcomm(fplog,cr,nstglobalcomm,ir);
1938 /*bGStatEveryStep = (nstglobalcomm == 1);*/
1939 bGStatEveryStep = FALSE;
1941 if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
1944 "To reduce the energy communication with nstlist = -1\n"
1945 "the neighbor list validity should not be checked at every step,\n"
1946 "this means that exact integration is not guaranteed.\n"
1947 "The neighbor list validity is checked after:\n"
1948 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
1949 "In most cases this will result in exact integration.\n"
1950 "This reduces the energy communication by a factor of 2 to 3.\n"
1951 "If you want less energy communication, set nstlist > 3.\n\n");
1954 if (bRerunMD || bFFscan)
1958 groups = &top_global->groups;
1960 /* Initial values */
1961 init_md(fplog,cr,ir,oenv,&t,&t0,&state_global->lambda,&lam0,
1962 nrnb,top_global,&upd,
1963 nfile,fnm,&outf,&mdebin,
1964 force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
1966 clear_mat(total_vir);
1968 /* Energy terms and groups */
1970 init_enerdata(top_global->groups.grps[egcENER].nr,ir->n_flambda,enerd);
1971 if (DOMAINDECOMP(cr))
1977 snew(f,top_global->natoms);
1980 /* Kinetic energy data */
1982 init_ekindata(fplog,top_global,&(ir->opts),ekind);
1983 /* needed for iteration of constraints */
1985 init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
1986 /* Copy the cos acceleration to the groups struct */
1987 ekind->cosacc.cos_accel = ir->cos_accel;
1989 gstat = global_stat_init(ir);
1992 /* Check for polarizable models and flexible constraints */
1993 shellfc = init_shell_flexcon(fplog,
1994 top_global,n_flexible_constraints(constr),
1995 (ir->bContinuation ||
1996 (DOMAINDECOMP(cr) && !MASTER(cr))) ?
1997 NULL : state_global->x);
2002 tMPI_Thread_mutex_lock(&deform_init_box_mutex);
2004 set_deform_reference_box(upd,
2005 deform_init_init_step_tpx,
2006 deform_init_box_tpx);
2008 tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
2013 double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
2014 if ((io > 2000) && MASTER(cr))
2016 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
2020 if (DOMAINDECOMP(cr)) {
2021 top = dd_init_local_top(top_global);
2024 dd_init_local_state(cr->dd,state_global,state);
2026 if (DDMASTER(cr->dd) && ir->nstfout) {
2027 snew(f_global,state_global->natoms);
2031 /* Initialize the particle decomposition and split the topology */
2032 top = split_system(fplog,top_global,ir,cr);
2034 pd_cg_range(cr,&fr->cg0,&fr->hcg);
2035 pd_at_range(cr,&a0,&a1);
2037 top = gmx_mtop_generate_local_top(top_global,ir);
2040 a1 = top_global->natoms;
2043 state = partdec_init_local_state(cr,state_global);
2046 atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
2049 set_vsite_top(vsite,top,mdatoms,cr);
2052 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols) {
2053 graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
2057 make_local_shells(cr,mdatoms,shellfc);
2060 if (ir->pull && PAR(cr)) {
2061 dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
2065 if (DOMAINDECOMP(cr))
2067 /* Distribute the charge groups over the nodes from the master node */
2068 dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
2069 state_global,top_global,ir,
2070 state,&f,mdatoms,top,fr,
2071 vsite,shellfc,constr,
2075 update_mdatoms(mdatoms,state->lambda);
2079 if (opt2bSet("-cpi",nfile,fnm))
2081 /* Update mdebin with energy history if appending to output files */
2082 if ( Flags & MD_APPENDFILES )
2084 restore_energyhistory_from_state(mdebin,&state_global->enerhist);
2088 /* We might have read an energy history from checkpoint,
2089 * free the allocated memory and reset the counts.
2091 done_energyhistory(&state_global->enerhist);
2092 init_energyhistory(&state_global->enerhist);
2095 /* Set the initial energy history in state by updating once */
2096 update_energyhistory(&state_global->enerhist,mdebin);
2099 if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG)) {
2100 /* Set the random state if we read a checkpoint file */
2101 set_stochd_state(upd,state);
2104 /* Initialize constraints */
2106 if (!DOMAINDECOMP(cr))
2107 set_constraints(constr,top,ir,mdatoms,cr);
2110 /* Check whether we have to GCT stuff */
2111 /* bTCR = ftp2bSet(efGCT,nfile,fnm);
2114 fprintf(stderr,"Will do General Coupling Theory!\n");
2116 gnx = top_global->mols.nr;
2118 for(i=0; (i<gnx); i++) {
2123 /* if (repl_ex_nst > 0 && MASTER(cr))
2124 repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
2125 repl_ex_nst,repl_ex_seed);*/
2127 if (!ir->bContinuation && !bRerunMD)
2129 if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
2131 /* Set the velocities of frozen particles to zero */
2132 for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
2134 for(m=0; m<DIM; m++)
2136 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
2146 /* Constrain the initial coordinates and velocities */
2147 do_constrain_first(fplog,constr,ir,mdatoms,state,f,
2148 graph,cr,nrnb,fr,top,shake_vir);
2152 /* Construct the virtual sites for the initial configuration */
2153 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
2154 top->idef.iparams,top->idef.il,
2155 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
2161 /* I'm assuming we need global communication the first time! MRS */
2162 cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
2163 | (bVV ? CGLO_PRESSURE:0)
2164 | (bVV ? CGLO_CONSTRAINT:0)
2165 | (bRerunMD ? CGLO_RERUNMD:0)
2166 | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
2168 bSumEkinhOld = FALSE;
2169 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2170 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2171 constr,NULL,FALSE,state->box,
2172 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
2173 if (ir->eI == eiVVAK) {
2174 /* a second call to get the half step temperature initialized as well */
2175 /* we do the same call as above, but turn the pressure off -- internally, this
2176 is recognized as a velocity verlet half-step kinetic energy calculation.
2177 This minimized excess variables, but perhaps loses some logic?*/
2179 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2180 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2181 constr,NULL,FALSE,state->box,
2182 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2183 cglo_flags &~ CGLO_PRESSURE);
2186 /* Calculate the initial half step temperature, and save the ekinh_old */
2187 if (!(Flags & MD_STARTFROMCPT))
2189 for(i=0; (i<ir->opts.ngtc); i++)
2191 copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
2196 enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
2197 and there is no previous step */
2199 temp0 = enerd->term[F_TEMP];
2201 /* if using an iterative algorithm, we need to create a working directory for the state. */
2204 bufstate = init_bufstate(state);
2208 snew(xcopy,state->natoms);
2209 snew(vcopy,state->natoms);
2210 copy_rvecn(state->x,xcopy,0,state->natoms);
2211 copy_rvecn(state->v,vcopy,0,state->natoms);
2212 copy_mat(state->box,boxcopy);
2215 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
2216 temperature control */
2217 trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
2221 if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
2224 "RMS relative constraint deviation after constraining: %.2e\n",
2225 constr_rmsd(constr,FALSE));
2227 fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
2230 fprintf(stderr,"starting md rerun '%s', reading coordinates from"
2231 " input trajectory '%s'\n\n",
2232 *(top_global->name),opt2fn("-rerun",nfile,fnm));
2235 fprintf(stderr,"Calculated time to finish depends on nsteps from "
2236 "run input file,\nwhich may not correspond to the time "
2237 "needed to process input trajectory.\n\n");
2243 fprintf(stderr,"starting mdrun '%s'\n",
2244 *(top_global->name));
2245 if (ir->nsteps >= 0)
2247 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
2251 sprintf(tbuf,"%s","infinite");
2253 if (ir->init_step > 0)
2255 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
2256 gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
2257 gmx_step_str(ir->init_step,sbuf2),
2258 ir->init_step*ir->delta_t);
2262 fprintf(stderr,"%s steps, %s ps.\n",
2263 gmx_step_str(ir->nsteps,sbuf),tbuf);
2266 fprintf(fplog,"\n");
2269 /* Set and write start time */
2270 runtime_start(runtime);
2271 print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
2272 wallcycle_start(wcycle,ewcRUN);
2274 fprintf(fplog,"\n");
2276 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
2277 /*#ifdef GMX_FAHCORE
2278 chkpt_ret=fcCheckPointParallel( cr->nodeid,
2280 if ( chkpt_ret == 0 )
2281 gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
2285 /***********************************************************
2287 * Loop over MD steps
2289 ************************************************************/
2291 /* if rerunMD then read coordinates and velocities from input trajectory */
2294 if (getenv("GMX_FORCE_UPDATE"))
2296 bForceUpdate = TRUE;
2299 bNotLastFrame = read_first_frame(oenv,&status,
2300 opt2fn("-rerun",nfile,fnm),
2301 &rerun_fr,TRX_NEED_X | TRX_READ_V);
2302 if (rerun_fr.natoms != top_global->natoms)
2305 "Number of atoms in trajectory (%d) does not match the "
2306 "run input file (%d)\n",
2307 rerun_fr.natoms,top_global->natoms);
2309 if (ir->ePBC != epbcNONE)
2313 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);
2315 if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
2317 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
2320 /* Set the shift vectors.
2321 * Necessary here when have a static box different from the tpr box.
2323 calc_shifts(rerun_fr.box,fr->shift_vec);
2327 /* loop over MD steps or if rerunMD to end of input trajectory */
2329 /* Skip the first Nose-Hoover integration when we get the state from tpx */
2330 bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
2331 bInitStep = bFirstStep && (bStateFromTPX || bVV);
2332 bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
2334 bSumEkinhOld = FALSE;
2337 init_global_signals(&gs,cr,ir,repl_ex_nst);
2339 step = ir->init_step;
2342 if (ir->nstlist == -1)
2344 init_nlistheuristics(&nlh,bGStatEveryStep,step);
2347 bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps));
2348 while (!bLastStep || (bRerunMD && bNotLastFrame)) {
2350 wallcycle_start(wcycle,ewcSTEP);
2352 GMX_MPE_LOG(ev_timestep1);
2355 if (rerun_fr.bStep) {
2356 step = rerun_fr.step;
2357 step_rel = step - ir->init_step;
2359 if (rerun_fr.bTime) {
2369 bLastStep = (step_rel == ir->nsteps);
2370 t = t0 + step*ir->delta_t;
2373 if (ir->efep != efepNO)
2375 if (bRerunMD && rerun_fr.bLambda && (ir->delta_lambda!=0))
2377 state_global->lambda = rerun_fr.lambda;
2381 state_global->lambda = lam0 + step*ir->delta_lambda;
2383 state->lambda = state_global->lambda;
2384 bDoDHDL = do_per_step(step,ir->nstdhdl);
2389 update_annealing_target_temp(&(ir->opts),t);
2394 if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
2396 for(i=0; i<state_global->natoms; i++)
2398 copy_rvec(rerun_fr.x[i],state_global->x[i]);
2402 for(i=0; i<state_global->natoms; i++)
2404 copy_rvec(rerun_fr.v[i],state_global->v[i]);
2409 for(i=0; i<state_global->natoms; i++)
2411 clear_rvec(state_global->v[i]);
2415 fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
2416 " Ekin, temperature and pressure are incorrect,\n"
2417 " the virial will be incorrect when constraints are present.\n"
2419 bRerunWarnNoV = FALSE;
2423 copy_mat(rerun_fr.box,state_global->box);
2424 copy_mat(state_global->box,state->box);
2426 if (vsite && (Flags & MD_RERUN_VSITE))
2428 if (DOMAINDECOMP(cr))
2430 gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
2434 /* Following is necessary because the graph may get out of sync
2435 * with the coordinates if we only have every N'th coordinate set
2437 mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
2438 shift_self(graph,state->box,state->x);
2440 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
2441 top->idef.iparams,top->idef.il,
2442 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
2445 unshift_self(graph,state->box,state->x);
2450 /* Stop Center of Mass motion */
2451 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
2453 /* Copy back starting coordinates in case we're doing a forcefield scan */
2456 for(ii=0; (ii<state->natoms); ii++)
2458 copy_rvec(xcopy[ii],state->x[ii]);
2459 copy_rvec(vcopy[ii],state->v[ii]);
2461 copy_mat(boxcopy,state->box);
2466 /* for rerun MD always do Neighbour Searching */
2467 bNS = (bFirstStep || ir->nstlist != 0);
2472 /* Determine whether or not to do Neighbour Searching and LR */
2473 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
2475 bNS = (bFirstStep || bExchanged || bNStList ||
2476 (ir->nstlist == -1 && nlh.nabnsb > 0));
2478 if (bNS && ir->nstlist == -1)
2480 set_nlistheuristics(&nlh,bFirstStep || bExchanged,step);
2484 /* < 0 means stop at next step, > 0 means stop at next NS step */
2485 if ( (gs.set[eglsSTOPCOND] < 0 ) ||
2486 ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
2491 /* Determine whether or not to update the Born radii if doing GB */
2492 bBornRadii=bFirstStep;
2493 if (ir->implicit_solvent && (step % ir->nstgbradii==0))
2498 do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
2499 do_verbose = bVerbose &&
2500 (step % stepout == 0 || bFirstStep || bLastStep);
2502 if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
2506 bMasterState = TRUE;
2510 bMasterState = FALSE;
2511 /* Correct the new box if it is too skewed */
2512 if (DYNAMIC_BOX(*ir))
2514 if (correct_box(fplog,step,state->box,graph))
2516 bMasterState = TRUE;
2519 if (DOMAINDECOMP(cr) && bMasterState)
2521 dd_collect_state(cr->dd,state,state_global);
2525 if (DOMAINDECOMP(cr))
2527 /* Repartition the domain decomposition */
2528 wallcycle_start(wcycle,ewcDOMDEC);
2529 dd_partition_system(fplog,step,cr,
2530 bMasterState,nstglobalcomm,
2531 state_global,top_global,ir,
2532 state,&f,mdatoms,top,fr,
2533 vsite,shellfc,constr,
2534 nrnb,wcycle,do_verbose);
2535 wallcycle_stop(wcycle,ewcDOMDEC);
2536 /* If using an iterative integrator, reallocate space to match the decomposition */
2540 if (MASTER(cr) && do_log && !bFFscan)
2542 print_ebin_header(fplog,step,t,state->lambda);
2545 if (ir->efep != efepNO)
2547 update_mdatoms(mdatoms,state->lambda);
2550 if (bRerunMD && rerun_fr.bV)
2553 /* We need the kinetic energy at minus the half step for determining
2554 * the full step kinetic energy and possibly for T-coupling.*/
2555 /* This may not be quite working correctly yet . . . . */
2556 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2557 wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
2558 constr,NULL,FALSE,state->box,
2559 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2560 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
2562 clear_mat(force_vir);
2564 /* Ionize the atoms if necessary */
2567 ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
2568 mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
2571 /* Update force field in ffscan program */
2574 if (update_forcefield(fplog,
2576 mdatoms->nr,state->x,state->box)) {
2577 if (gmx_parallel_env_initialized())
2585 GMX_MPE_LOG(ev_timestep2);
2587 /* We write a checkpoint at this MD step when:
2588 * either at an NS step when we signalled through gs,
2589 * or at the last step (but not when we do not want confout),
2590 * but never at the first step or with rerun.
2592 /* bCPT = (((gs.set[eglsCHKPT] && bNS) ||
2593 (bLastStep && (Flags & MD_CONFOUT))) &&
2594 step > ir->init_step && !bRerunMD);
2597 gs.set[eglsCHKPT] = 0;
2600 /* Determine the energy and pressure:
2601 * at nstcalcenergy steps and at energy output steps (set below).
2603 bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
2604 bCalcEnerPres = bNstEner;
2606 /* Do we need global communication ? */
2607 bGStat = (bCalcEnerPres || bStopCM ||
2608 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
2610 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
2612 if (do_ene || do_log)
2614 bCalcEnerPres = TRUE;
2618 /* these CGLO_ options remain the same throughout the iteration */
2619 cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
2620 (bStopCM ? CGLO_STOPCM : 0) |
2621 (bGStat ? CGLO_GSTAT : 0)
2624 force_flags = (GMX_FORCE_STATECHANGED |
2625 ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
2626 GMX_FORCE_ALLFORCES |
2627 (bNStList ? GMX_FORCE_DOLR : 0) |
2629 (bCalcEnerPres ? GMX_FORCE_VIRIAL : 0) |
2630 (bDoDHDL ? GMX_FORCE_DHDL : 0)
2635 /* Now is the time to relax the shells */
2636 count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
2638 bStopCM,top,top_global,
2640 state,f,force_vir,mdatoms,
2641 nrnb,wcycle,graph,groups,
2642 shellfc,fr,bBornRadii,t,mu_tot,
2643 state->natoms,&bConverged,vsite,
2654 /* The coordinates (x) are shifted (to get whole molecules)
2656 * This is parallellized as well, and does communication too.
2657 * Check comments in sim_util.c
2660 do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
2661 state->box,state->x,&state->hist,
2662 f,force_vir,mdatoms,enerd,fcd,
2663 state->lambda,graph,
2664 fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
2665 (bNS ? GMX_FORCE_NS : 0) | force_flags);
2668 GMX_BARRIER(cr->mpi_comm_mygroup);
2672 mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
2673 mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
2676 if (bTCR && bFirstStep)
2678 tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
2679 fprintf(fplog,"Done init_coupling\n");
2683 /* ############### START FIRST UPDATE HALF-STEP ############### */
2685 if (bVV && !bStartingFromCpt && !bRerunMD)
2691 /* if using velocity verlet with full time step Ekin,
2692 * take the first half step only to compute the
2693 * virial for the first step. From there,
2694 * revert back to the initial coordinates
2695 * so that the input is actually the initial step.
2697 copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
2700 /* this is for NHC in the Ekin(t+dt/2) version of vv */
2703 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
2706 if (ir->eI == eiVVAK)
2708 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
2711 update_coords(fplog,step,ir,mdatoms,state,
2712 f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2713 ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
2714 cr,nrnb,constr,&top->idef);
2718 gmx_iterate_init(&iterate,bIterations && !bInitStep);
2720 /* for iterations, we save these vectors, as we will be self-consistently iterating
2722 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
2724 /* save the state */
2725 if (bIterations && iterate.bIterate) {
2726 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
2730 bFirstIterate = TRUE;
2731 while (bFirstIterate || (bIterations && iterate.bIterate))
2733 if (bIterations && iterate.bIterate)
2735 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
2736 if (bFirstIterate && bTrotter)
2738 /* The first time through, we need a decent first estimate
2739 of veta(t+dt) to compute the constraints. Do
2740 this by computing the box volume part of the
2741 trotter integration at this time. Nothing else
2742 should be changed by this routine here. If
2743 !(first time), we start with the previous value
2746 veta_save = state->veta;
2747 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
2748 vetanew = state->veta;
2749 state->veta = veta_save;
2754 if ( !bRerunMD || rerun_fr.bV || bForceUpdate) { /* Why is rerun_fr.bV here? Unclear. */
2757 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2758 &top->idef,shake_vir,NULL,
2759 cr,nrnb,wcycle,upd,constr,
2760 bInitStep,TRUE,bCalcEnerPres,vetanew);
2762 if (!bOK && !bFFscan)
2764 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2769 { /* Need to unshift here if a do_force has been
2770 called in the previous step */
2771 unshift_self(graph,state->box,state->x);
2776 /* if VV, compute the pressure and constraints */
2777 /* if VV2, the pressure and constraints only if using pressure control.*/
2778 bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir));
2779 bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));
2780 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2781 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2782 constr,NULL,FALSE,state->box,
2783 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2786 | (bTemp ? CGLO_TEMPERATURE:0)
2787 | (bPres ? CGLO_PRESSURE : 0)
2788 | (bPres ? CGLO_CONSTRAINT : 0)
2789 | (iterate.bIterate ? CGLO_ITERATE : 0)
2790 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
2794 /* explanation of above:
2795 a) We compute Ekin at the full time step
2796 if 1) we are using the AveVel Ekin, and it's not the
2797 initial step, or 2) if we are using AveEkin, but need the full
2798 time step kinetic energy for the pressure.
2799 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
2800 EkinAveVel because it's needed for the pressure */
2802 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
2803 if (bVV && !bInitStep)
2805 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ2);
2809 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
2810 state->veta,&vetanew))
2814 bFirstIterate = FALSE;
2817 if (bTrotter && !bInitStep) {
2818 copy_mat(shake_vir,state->svir_prev);
2819 copy_mat(force_vir,state->fvir_prev);
2820 if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
2821 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
2822 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
2823 enerd->term[F_EKIN] = trace(ekind->ekin);
2826 /* if it's the initial step, we performed this first step just to get the constraint virial */
2827 if (bInitStep && ir->eI==eiVV) {
2828 copy_rvecn(cbuf,state->v,0,state->natoms);
2831 if (fr->bSepDVDL && fplog && do_log)
2833 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
2835 enerd->term[F_DHDL_CON] += dvdl;
2837 GMX_MPE_LOG(ev_timestep1);
2841 /* MRS -- now done iterating -- compute the conserved quantity */
2844 if (IR_NVT_TROTTER(ir) || IR_NPT_TROTTER(ir))
2847 NPT_energy(ir,state,&MassQ);
2848 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
2850 last_conserved -= enerd->term[F_DISPCORR];
2854 last_ekin = enerd->term[F_EKIN]; /* does this get preserved through checkpointing? */
2858 /* ######## END FIRST UPDATE STEP ############## */
2859 /* ######## If doing VV, we now have v(dt) ###### */
2861 /* ################## START TRAJECTORY OUTPUT ################# */
2863 /* Now we have the energies and forces corresponding to the
2864 * coordinates at time t. We must output all of this before
2866 * for RerunMD t is read from input trajectory
2868 GMX_MPE_LOG(ev_output_start);
2871 if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
2872 if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
2873 if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
2874 if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
2875 /* if (bCPT) { mdof_flags |= MDOF_CPT; };*/
2879 fcReportProgress( ir->nsteps, step );
2883 /* Enforce writing positions and velocities at end of run */
2884 mdof_flags |= (MDOF_X | MDOF_V);
2886 /* sync bCPT and fc record-keeping */
2887 /* if (bCPT && MASTER(cr))
2888 fcRequestCheckPoint();*/
2891 if (mdof_flags != 0)
2893 wallcycle_start(wcycle,ewcTRAJ);
2896 if (state->flags & (1<<estLD_RNG))
2898 get_stochd_state(upd,state);
2904 state_global->ekinstate.bUpToDate = FALSE;
2908 update_ekinstate(&state_global->ekinstate,ekind);
2909 state_global->ekinstate.bUpToDate = TRUE;
2911 update_energyhistory(&state_global->enerhist,mdebin);
2914 write_traj(fplog,cr,outf,mdof_flags,top_global,
2915 step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
2922 if (bLastStep && step_rel == ir->nsteps &&
2923 (Flags & MD_CONFOUT) && MASTER(cr) &&
2924 !bRerunMD && !bFFscan)
2926 /* x and v have been collected in write_traj,
2927 * because a checkpoint file will always be written
2930 fprintf(stderr,"\nWriting final coordinates.\n");
2931 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
2934 /* Make molecules whole only for confout writing */
2935 do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
2937 /* write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
2938 *top_global->name,top_global,
2939 state_global->x,state_global->v,
2940 ir->ePBC,state->box);*/
2943 wallcycle_stop(wcycle,ewcTRAJ);
2945 GMX_MPE_LOG(ev_output_finish);
2947 /* kludge -- virial is lost with restart for NPT control. Must restart */
2948 if (bStartingFromCpt && bVV)
2950 copy_mat(state->svir_prev,shake_vir);
2951 copy_mat(state->fvir_prev,force_vir);
2953 /* ################## END TRAJECTORY OUTPUT ################ */
2955 /* Determine the pressure:
2956 * always when we want exact averages in the energy file,
2957 * at ns steps when we have pressure coupling,
2958 * otherwise only at energy output steps (set below).
2961 bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
2962 bCalcEnerPres = bNstEner;
2964 /* Do we need global communication ? */
2965 bGStat = (bGStatEveryStep || bStopCM || bNS ||
2966 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
2968 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
2970 if (do_ene || do_log)
2972 bCalcEnerPres = TRUE;
2976 /* Determine the wallclock run time up till now */
2977 run_time = gmx_gettime() - (double)runtime->real;
2979 /* Check whether everything is still allright */
2980 if (((int)gmx_get_stop_condition() > handled_stop_condition)
2986 /* this is just make gs.sig compatible with the hack
2987 of sending signals around by MPI_Reduce with together with
2989 if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
2990 gs.sig[eglsSTOPCOND]=1;
2991 if ( gmx_get_stop_condition() == gmx_stop_cond_next )
2992 gs.sig[eglsSTOPCOND]=-1;
2993 /* < 0 means stop at next step, > 0 means stop at next NS step */
2997 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2998 gmx_get_signal_name(),
2999 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
3003 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
3004 gmx_get_signal_name(),
3005 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
3007 handled_stop_condition=(int)gmx_get_stop_condition();
3009 else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
3010 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
3011 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
3013 /* Signal to terminate the run */
3014 gs.sig[eglsSTOPCOND] = 1;
3017 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
3019 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
3022 if (bResetCountersHalfMaxH && MASTER(cr) &&
3023 run_time > max_hours*60.0*60.0*0.495)
3025 gs.sig[eglsRESETCOUNTERS] = 1;
3028 if (ir->nstlist == -1 && !bRerunMD)
3030 /* When bGStatEveryStep=FALSE, global_stat is only called
3031 * when we check the atom displacements, not at NS steps.
3032 * This means that also the bonded interaction count check is not
3033 * performed immediately after NS. Therefore a few MD steps could
3034 * be performed with missing interactions.
3035 * But wrong energies are never written to file,
3036 * since energies are only written after global_stat
3039 if (step >= nlh.step_nscheck)
3041 nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
3042 nlh.scale_tot,state->x);
3046 /* This is not necessarily true,
3047 * but step_nscheck is determined quite conservatively.
3053 /* In parallel we only have to check for checkpointing in steps
3054 * where we do global communication,
3055 * otherwise the other nodes don't know.
3057 if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
3060 run_time >= nchkpt*cpt_period*60.0)) &&
3061 gs.set[eglsCHKPT] == 0)
3063 gs.sig[eglsCHKPT] = 1;
3068 gmx_iterate_init(&iterate,bIterations);
3071 /* for iterations, we save these vectors, as we will be redoing the calculations */
3072 if (bIterations && iterate.bIterate)
3074 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
3076 bFirstIterate = TRUE;
3077 while (bFirstIterate || (bIterations && iterate.bIterate))
3079 /* We now restore these vectors to redo the calculation with improved extended variables */
3082 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
3085 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
3086 so scroll down for that logic */
3088 /* ######### START SECOND UPDATE STEP ################# */
3089 GMX_MPE_LOG(ev_update_start);
3091 if (!bRerunMD || rerun_fr.bV || bForceUpdate)
3093 wallcycle_start(wcycle,ewcUPDATE);
3095 /* Box is changed in update() when we do pressure coupling,
3096 * but we should still use the old box for energy corrections and when
3097 * writing it to the energy file, so it matches the trajectory files for
3098 * the same timestep above. Make a copy in a separate array.
3100 copy_mat(state->box,lastbox);
3101 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
3104 if (bIterations && iterate.bIterate)
3112 /* we use a new value of scalevir to converge the iterations faster */
3113 scalevir = tracevir/trace(shake_vir);
3115 msmul(shake_vir,scalevir,shake_vir);
3116 m_add(force_vir,shake_vir,total_vir);
3117 clear_mat(shake_vir);
3119 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ3);
3121 /* We can only do Berendsen coupling after we have summed
3122 * the kinetic energy or virial. Since the happens
3123 * in global_state after update, we should only do it at
3124 * step % nstlist = 1 with bGStatEveryStep=FALSE.
3127 if (ir->eI != eiVVAK)
3129 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
3131 update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
3136 /* velocity half-step update */
3137 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
3138 ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,cr,nrnb,constr,&top->idef);
3141 /* Above, initialize just copies ekinh into ekin,
3142 * it doesn't copy position (for VV),
3143 * and entire integrator for MD.
3148 copy_rvecn(state->x,cbuf,0,state->natoms);
3151 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
3152 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
3153 wallcycle_stop(wcycle,ewcUPDATE);
3155 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
3156 &top->idef,shake_vir,force_vir,
3157 cr,nrnb,wcycle,upd,constr,
3158 bInitStep,FALSE,bCalcEnerPres,state->veta);
3162 /* erase F_EKIN and F_TEMP here? */
3163 /* just compute the kinetic energy at the half step to perform a trotter step */
3164 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
3165 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
3166 constr,NULL,FALSE,lastbox,
3167 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
3168 cglo_flags | CGLO_TEMPERATURE | CGLO_CONSTRAINT
3170 wallcycle_start(wcycle,ewcUPDATE);
3171 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ4);
3172 /* now we know the scaling, we can compute the positions again again */
3173 copy_rvecn(cbuf,state->x,0,state->natoms);
3175 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
3176 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
3177 wallcycle_stop(wcycle,ewcUPDATE);
3179 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
3180 /* are the small terms in the shake_vir here due
3181 * to numerical errors, or are they important
3182 * physically? I'm thinking they are just errors, but not completely sure.
3183 * For now, will call without actually constraining, constr=NULL*/
3184 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
3185 &top->idef,tmp_vir,force_vir,
3186 cr,nrnb,wcycle,upd,NULL,
3187 bInitStep,FALSE,bCalcEnerPres,state->veta);
3189 if (!bOK && !bFFscan)
3191 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
3194 if (fr->bSepDVDL && fplog && do_log)
3196 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
3198 enerd->term[F_DHDL_CON] += dvdl;
3202 /* Need to unshift here */
3203 unshift_self(graph,state->box,state->x);
3206 GMX_BARRIER(cr->mpi_comm_mygroup);
3207 GMX_MPE_LOG(ev_update_finish);
3211 wallcycle_start(wcycle,ewcVSITECONSTR);
3214 shift_self(graph,state->box,state->x);
3216 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
3217 top->idef.iparams,top->idef.il,
3218 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
3222 unshift_self(graph,state->box,state->x);
3224 wallcycle_stop(wcycle,ewcVSITECONSTR);
3227 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
3228 if (ir->nstlist == -1 && bFirstIterate)
3230 gs.sig[eglsNABNSB] = nlh.nabnsb;
3232 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
3233 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
3235 bFirstIterate ? &gs : NULL,(step % gs.nstms == 0),
3237 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
3239 | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
3240 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
3241 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
3242 | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
3243 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
3246 if (ir->nstlist == -1 && bFirstIterate)
3248 nlh.nabnsb = gs.set[eglsNABNSB];
3249 gs.set[eglsNABNSB] = 0;
3251 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
3252 /* ############# END CALC EKIN AND PRESSURE ################# */
3254 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
3255 the virial that should probably be addressed eventually. state->veta has better properies,
3256 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
3257 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
3260 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
3261 trace(shake_vir),&tracevir))
3265 bFirstIterate = FALSE;
3268 update_box(fplog,step,ir,mdatoms,state,graph,f,
3269 ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
3271 /* ################# END UPDATE STEP 2 ################# */
3272 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
3274 /* The coordinates (x) were unshifted in update */
3275 /* if (bFFscan && (shellfc==NULL || bConverged))
3277 if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
3279 &(top_global->mols),mdatoms->massT,pres))
3281 if (gmx_parallel_env_initialized())
3285 fprintf(stderr,"\n");
3291 /* We will not sum ekinh_old,
3292 * so signal that we still have to do it.
3294 bSumEkinhOld = TRUE;
3299 /* Only do GCT when the relaxation of shells (minimization) has converged,
3300 * otherwise we might be coupling to bogus energies.
3301 * In parallel we must always do this, because the other sims might
3305 /* Since this is called with the new coordinates state->x, I assume
3306 * we want the new box state->box too. / EL 20040121
3308 /* do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
3310 mdatoms,&(top->idef),mu_aver,
3311 top_global->mols.nr,cr,
3312 state->box,total_vir,pres,
3313 mu_tot,state->x,f,bConverged);
3317 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
3319 sum_dhdl(enerd,state->lambda,ir);
3320 /* use the directly determined last velocity, not actually the averaged half steps */
3321 if (bTrotter && ir->eI==eiVV)
3323 enerd->term[F_EKIN] = last_ekin;
3325 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
3334 if (IR_NVT_TROTTER(ir)) {
3335 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + last_conserved;
3337 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] +
3338 NPT_energy(ir,state,&MassQ);
3342 enerd->term[F_ECONSERVED] =
3343 enerd->term[F_ETOT] + vrescale_energy(&(ir->opts),
3344 state->therm_integral);
3350 /* Check for excessively large energies */
3354 real etot_max = 1e200;
3356 real etot_max = 1e30;
3358 if (fabs(enerd->term[F_ETOT]) > etot_max)
3360 fprintf(stderr,"Energy too large (%g), giving up\n",
3361 enerd->term[F_ETOT]);
3364 /* ######### END PREPARING EDR OUTPUT ########### */
3366 /* Time for performance */
3367 if (((step % stepout) == 0) || bLastStep)
3369 runtime_upd_proc(runtime);
3375 gmx_bool do_dr,do_or;
3377 if (!(bStartingFromCpt && (EI_VV(ir->eI))))
3381 upd_mdebin(mdebin,bDoDHDL,TRUE,
3382 t,mdatoms->tmass,enerd,state,lastbox,
3383 shake_vir,force_vir,total_vir,pres,
3384 ekind,mu_tot,constr);
3388 upd_mdebin_step(mdebin);
3391 do_dr = do_per_step(step,ir->nstdisreout);
3392 do_or = do_per_step(step,ir->nstorireout);
3394 print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
3396 eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
3398 if (ir->ePull != epullNO)
3400 pull_print_output(ir->pull,step,t);
3403 if (do_per_step(step,ir->nstlog))
3405 if(fflush(fplog) != 0)
3407 gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of quota?");
3413 /* Remaining runtime */
3414 if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
3418 fprintf(stderr,"\n");
3420 print_time(stderr,runtime,step,ir,cr);
3423 /* Set new positions for the group to embed */
3429 } else if (step_rel<=(it_xy+it_z))
3433 resize(ins_at,r_ins,state_global->x,pos_ins,fac);
3436 /* Replica exchange */
3437 /* bExchanged = FALSE;
3438 if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
3439 do_per_step(step,repl_ex_nst))
3441 bExchanged = replica_exchange(fplog,cr,repl_ex,
3442 state_global,enerd->term,
3445 if (bExchanged && PAR(cr))
3447 if (DOMAINDECOMP(cr))
3449 dd_partition_system(fplog,step,cr,TRUE,1,
3450 state_global,top_global,ir,
3451 state,&f,mdatoms,top,fr,
3452 vsite,shellfc,constr,
3457 bcast_state(cr,state,FALSE);
3463 bStartingFromCpt = FALSE;
3465 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
3466 /* With all integrators, except VV, we need to retain the pressure
3467 * at the current step for coupling at the next step.
3469 if ((state->flags & (1<<estPRES_PREV)) &&
3471 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
3473 /* Store the pressure in t_state for pressure coupling
3474 * at the next MD step.
3476 copy_mat(pres,state->pres_prev);
3479 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
3483 /* read next frame from input trajectory */
3484 bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
3487 if (!bRerunMD || !rerun_fr.bStep)
3489 /* increase the MD step number */
3494 cycles = wallcycle_stop(wcycle,ewcSTEP);
3495 if (DOMAINDECOMP(cr) && wcycle)
3497 dd_cycles_add(cr->dd,cycles,ddCyclStep);
3500 if (step_rel == wcycle_get_reset_counters(wcycle) ||
3501 gs.set[eglsRESETCOUNTERS] != 0)
3503 /* Reset all the counters related to performance over the run */
3504 reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime);
3505 wcycle_set_reset_counters(wcycle,-1);
3506 bResetCountersHalfMaxH = FALSE;
3507 gs.set[eglsRESETCOUNTERS] = 0;
3510 /* End of main MD loop */
3512 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
3513 *top_global->name,top_global,
3514 state_global->x,state_global->v,
3515 ir->ePBC,state->box);
3518 runtime_end(runtime);
3525 if (!(cr->duty & DUTY_PME))
3527 /* Tell the PME only node to finish */
3533 if (ir->nstcalcenergy > 0 && !bRerunMD)
3535 print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
3536 eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
3544 if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
3546 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)));
3547 fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
3550 if (shellfc && fplog)
3552 fprintf(fplog,"Fraction of iterations that converged: %.2f %%\n",
3553 (nconverged*100.0)/step_rel);
3554 fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
3558 /* if (repl_ex_nst > 0 && MASTER(cr))
3560 print_replica_exchange_statistics(fplog,repl_ex);
3563 runtime->nsteps_done = step_rel;
3569 int mdrunner_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
3570 const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
3572 ivec ddxyz,int dd_node_order,real rdd,real rconstr,
3573 const char *dddlb_opt,real dlb_scale,
3574 const char *ddcsx,const char *ddcsy,const char *ddcsz,
3575 int nstepout,int resetstep,int nmultisim,int repl_ex_nst,int repl_ex_seed,
3576 real pforce,real cpt_period,real max_hours,
3577 const char *deviceOptions,
3578 unsigned long Flags,
3579 real xy_fac, real xy_max, real z_fac, real z_max,
3580 int it_xy, int it_z, real probe_rad, int low_up_rm,
3581 int pieces, gmx_bool bALLOW_ASYMMETRY, int maxwarn)
3583 double nodetime=0,realtime;
3584 t_inputrec *inputrec;
3585 t_state *state=NULL;
3588 int npme_major,npme_minor;
3591 gmx_mtop_t *mtop=NULL;
3592 t_mdatoms *mdatoms=NULL;
3593 t_forcerec *fr=NULL;
3596 gmx_pme_t *pmedata=NULL;
3597 gmx_vsite_t *vsite=NULL;
3598 gmx_constr_t constr;
3599 int i,m,nChargePerturbed=-1,status,nalloc;
3601 gmx_wallcycle_t wcycle;
3602 gmx_bool bReadRNG,bReadEkin;
3604 gmx_runtime_t runtime;
3606 gmx_large_int_t reset_counters;
3607 gmx_edsam_t ed=NULL;
3608 t_commrec *cr_old=cr;
3609 int nthreads=1,nthreads_requested=1;
3613 int rm_bonded_at,fr_id,fr_i=0,tmp_id,warn=0;
3614 int ng,j,max_lip_rm,ins_grp_id,ins_nat,mem_nat,ntype,lip_rm,tpr_version;
3615 real xy_step=0,z_step=0;
3617 rvec *r_ins=NULL,fac;
3618 t_block *ins_at,*rest_at;
3622 gmx_groups_t *groups;
3623 gmx_bool bExcl=FALSE;
3626 char **piecename=NULL;
3628 /* CAUTION: threads may be started later on in this function, so
3629 cr doesn't reflect the final parallel state right now */
3633 if (bVerbose && SIMMASTER(cr))
3635 fprintf(stderr,"Getting Loaded...\n");
3638 if (Flags & MD_APPENDFILES)
3646 /* Read (nearly) all data required for the simulation */
3647 read_tpx_state(ftp2fn(efTPX,nfile,fnm),inputrec,state,NULL,mtop);
3649 /* NOW the threads will be started: */
3653 /* END OF CAUTION: cr is now reliable */
3657 /* now broadcast everything to the non-master nodes/threads: */
3658 init_parallel(fplog, cr, inputrec, mtop);
3660 /* now make sure the state is initialized and propagated */
3661 set_state_entries(state,inputrec,cr->nnodes);
3663 if (can_use_allvsall(inputrec,mtop,TRUE,cr,fplog))
3665 /* All-vs-all loops do not work with domain decomposition */
3666 Flags |= MD_PARTDEC;
3669 if (!EEL_PME(inputrec->coulombtype) || (Flags & MD_PARTDEC))
3678 tpr_version = get_tpr_version(ftp2fn(efTPX,nfile,fnm));
3680 gmx_fatal(FARGS,"Version of *.tpr file to old (%d). Rerun grompp with gromacs VERSION 4.0.3 or newer.\n",tpr_version);
3682 if( inputrec->eI != eiMD )
3683 gmx_input("Change integrator to md in mdp file.");
3686 gmx_input("Sorry, parallel g_membed is not yet fully functrional.");
3688 groups=&(mtop->groups);
3690 atoms=gmx_mtop_global_atoms(mtop);
3692 fprintf(stderr,"\nSelect a group to embed in the membrane:\n");
3693 get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),1,&(ins_at->nr),&(ins_at->index),&ins);
3694 ins_grp_id = search_string(ins,groups->ngrpname,(groups->grpname));
3695 fprintf(stderr,"\nSelect a group to embed %s into (e.g. the membrane):\n",ins);
3696 get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),1,&(mem_p->mem_at.nr),&(mem_p->mem_at.index),&(mem_p->name));
3698 pos_ins->pieces=pieces;
3699 snew(pos_ins->nidx,pieces);
3700 snew(pos_ins->subindex,pieces);
3701 snew(piecename,pieces);
3704 fprintf(stderr,"\nSelect pieces to embed:\n");
3705 get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),pieces,pos_ins->nidx,pos_ins->subindex,piecename);
3709 /*use whole embedded group*/
3710 snew(pos_ins->nidx,1);
3711 snew(pos_ins->subindex,1);
3712 pos_ins->nidx[0]=ins_at->nr;
3713 pos_ins->subindex[0]=ins_at->index;
3716 if(probe_rad<0.2199999)
3719 fprintf(stderr,"\nWarning %d:\nA probe radius (-rad) smaller than 0.2 can result in overlap between waters "
3720 "and the group to embed, which will result in Lincs errors etc.\nIf you are sure, you can increase maxwarn.\n\n",warn);
3723 if(xy_fac<0.09999999)
3726 fprintf(stderr,"\nWarning %d:\nThe initial size of %s is probably too smal.\n"
3727 "If you are sure, you can increase maxwarn.\n\n",warn,ins);
3733 fprintf(stderr,"\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d) is probably too small.\n"
3734 "Increase -nxy or, if you are sure, you can increase maxwarn.\n\n",warn,ins,it_xy);
3737 if( (it_z<100) && ( z_fac<0.99999999 || z_fac>1.0000001) )
3740 fprintf(stderr,"\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d) is probably too small.\n"
3741 "Increase -nz or, if you are sure, you can increase maxwarn.\n\n",warn,ins,it_z);
3744 if(it_xy+it_z>inputrec->nsteps)
3747 fprintf(stderr,"\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the number of steps in the tpr.\n"
3748 "If you are sure, you can increase maxwarn.\n\n",warn);
3752 if( inputrec->opts.ngfrz==1)
3753 gmx_fatal(FARGS,"You did not specify \"%s\" as a freezegroup.",ins);
3754 for(i=0;i<inputrec->opts.ngfrz;i++)
3756 tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i];
3757 if(ins_grp_id==tmp_id)
3764 gmx_fatal(FARGS,"\"%s\" not as freezegroup defined in the mdp-file.",ins);
3767 if( inputrec->opts.nFreeze[fr_i][i] != 1)
3768 gmx_fatal(FARGS,"freeze dimensions for %s are not Y Y Y\n",ins);
3770 ng = groups->grps[egcENER].nr;
3772 gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
3778 if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL)
3781 if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) || (groups->grps[egcENER].nm_ind[j] != ins_grp_id) )
3782 gmx_fatal(FARGS,"Energy exclusions \"%s\" and \"%s\" do not match the group to embed \"%s\"",
3783 *groups->grpname[groups->grps[egcENER].nm_ind[i]],
3784 *groups->grpname[groups->grps[egcENER].nm_ind[j]],ins);
3789 gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in the freeze group");
3791 /* Set all atoms in box*/
3792 /*set_inbox(state->natoms,state->x);*/
3794 /* Guess the area the protein will occupy in the membrane plane Calculate area per lipid*/
3796 ins_nat = init_ins_at(ins_at,rest_at,state,pos_ins,groups,ins_grp_id,xy_max);
3797 /* Check moleculetypes in insertion group */
3798 check_types(ins_at,rest_at,mtop);
3800 mem_nat = init_mem_at(mem_p,mtop,state->x,state->box,pos_ins);
3802 prot_area = est_prot_area(pos_ins,state->x,ins_at,mem_p);
3803 if ( (prot_area>7.5) && ( (state->box[XX][XX]*state->box[YY][YY]-state->box[XX][YY]*state->box[YY][XX])<50) )
3806 fprintf(stderr,"\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
3807 "This might cause pressure problems during the growth phase. Just try with\n"
3808 "current setup (-maxwarn + 1), but if pressure problems occur, lower the\n"
3809 "compressibility in the mdp-file or use no pressure coupling at all.\n\n",warn);
3812 gmx_fatal(FARGS,"Too many warnings.\n");
3814 printf("The estimated area of the protein in the membrane is %.3f nm^2\n",prot_area);
3815 printf("\nThere are %d lipids in the membrane part that overlaps the protein.\nThe area per lipid is %.4f nm^2.\n",mem_p->nmol,mem_p->lip_area);
3817 /* Maximum number of lipids to be removed*/
3818 max_lip_rm=(int)(2*prot_area/mem_p->lip_area);
3819 printf("Maximum number of lipids that will be removed is %d.\n",max_lip_rm);
3821 printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
3822 "This resizing will be done with respect to the geometrical center of all protein atoms\n"
3823 "that span the membrane region, i.e. z between %.3f and %.3f\n\n",xy_fac,z_fac,mem_p->zmin,mem_p->zmax);
3825 /* resize the protein by xy and by z if necessary*/
3826 snew(r_ins,ins_at->nr);
3827 init_resize(ins_at,r_ins,pos_ins,mem_p,state->x,bALLOW_ASYMMETRY);
3828 fac[0]=fac[1]=xy_fac;
3831 xy_step =(xy_max-xy_fac)/(double)(it_xy);
3832 z_step =(z_max-z_fac)/(double)(it_z-1);
3834 resize(ins_at,r_ins,state->x,pos_ins,fac);
3836 /* remove overlapping lipids and water from the membrane box*/
3837 /*mark molecules to be removed*/
3839 set_pbc(pbc,inputrec->ePBC,state->box);
3842 lip_rm = gen_rm_list(rm_p,ins_at,rest_at,pbc,mtop,state->x, r_ins, mem_p,pos_ins,probe_rad,low_up_rm,bALLOW_ASYMMETRY);
3843 lip_rm -= low_up_rm;
3846 for(i=0;i<rm_p->nr;i++)
3847 fprintf(fplog,"rm mol %d\n",rm_p->mol[i]);
3849 for(i=0;i<mtop->nmolblock;i++)
3852 for(j=0;j<rm_p->nr;j++)
3853 if(rm_p->block[j]==i)
3855 printf("Will remove %d %s molecules\n",ntype,*(mtop->moltype[mtop->molblock[i].type].name));
3858 if(lip_rm>max_lip_rm)
3861 fprintf(stderr,"\nWarning %d:\nTrying to remove a larger lipid area than the estimated protein area\n"
3862 "Try making the -xyinit resize factor smaller. If you are sure about this increase maxwarn.\n\n",warn);
3865 /*remove all lipids and waters overlapping and update all important structures*/
3866 rm_group(inputrec,groups,mtop,rm_p,state,ins_at,pos_ins);
3868 rm_bonded_at = rm_bonded(ins_at,mtop);
3869 if (rm_bonded_at != ins_at->nr)
3871 fprintf(stderr,"Warning: The number of atoms for which the bonded interactions are removed is %d, "
3872 "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
3873 "molecule type as atoms that are not to be embedded.\n",rm_bonded_at,ins_at->nr);
3877 gmx_fatal(FARGS,"Too many warnings.\nIf you are sure these warnings are harmless, you can increase -maxwarn");
3881 if (ftp2bSet(efTOP,nfile,fnm))
3882 top_update(opt2fn("-p",nfile,fnm),ins,rm_p,mtop);
3890 fcRegisterSteps(inputrec->nsteps,inputrec->init_step);
3893 /* NMR restraints must be initialized before load_checkpoint,
3894 * since with time averaging the history is added to t_state.
3895 * For proper consistency check we therefore need to extend
3897 * So the PME-only nodes (if present) will also initialize
3898 * the distance restraints.
3902 /* This needs to be called before read_checkpoint to extend the state */
3903 init_disres(fplog,mtop,inputrec,cr,Flags & MD_PARTDEC,fcd,state);
3905 if (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0)
3907 if (PAR(cr) && !(Flags & MD_PARTDEC))
3909 gmx_fatal(FARGS,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
3911 /* Orientation restraints */
3914 init_orires(fplog,mtop,state->x,inputrec,cr->ms,&(fcd->orires),
3919 if (DEFORM(*inputrec))
3921 /* Store the deform reference box before reading the checkpoint */
3924 copy_mat(state->box,box);
3928 gmx_bcast(sizeof(box),box,cr);
3930 /* Because we do not have the update struct available yet
3931 * in which the reference values should be stored,
3932 * we store them temporarily in static variables.
3933 * This should be thread safe, since they are only written once
3934 * and with identical values.
3936 /* deform_init_init_step_tpx = inputrec->init_step;*/
3937 /* copy_mat(box,deform_init_box_tpx);*/
3940 if (opt2bSet("-cpi",nfile,fnm))
3942 /* Check if checkpoint file exists before doing continuation.
3943 * This way we can use identical input options for the first and subsequent runs...
3945 if( gmx_fexist_master(opt2fn_master("-cpi",nfile,fnm,cr),cr) )
3947 load_checkpoint(opt2fn_master("-cpi",nfile,fnm,cr),&fplog,
3948 cr,Flags & MD_PARTDEC,ddxyz,
3949 inputrec,state,&bReadRNG,&bReadEkin,
3950 (Flags & MD_APPENDFILES));
3954 Flags |= MD_READ_RNG;
3958 Flags |= MD_READ_EKIN;
3963 if ((MASTER(cr) || (Flags & MD_SEPPOT)) && (Flags & MD_APPENDFILES))
3965 gmx_log_open(ftp2fn(efLOG,nfile,fnm),cr,!(Flags & MD_SEPPOT),
3971 copy_mat(state->box,box);
3976 gmx_bcast(sizeof(box),box,cr);
3979 if (bVerbose && SIMMASTER(cr))
3981 fprintf(stderr,"Loaded with Money\n\n");
3984 if (PAR(cr) && !((Flags & MD_PARTDEC) || EI_TPI(inputrec->eI)))
3986 cr->dd = init_domain_decomposition(fplog,cr,Flags,ddxyz,rdd,rconstr,
3987 dddlb_opt,dlb_scale,
3991 &ddbox,&npme_major,&npme_minor);
3993 make_dd_communicators(fplog,cr,dd_node_order);
3995 /* Set overallocation to avoid frequent reallocation of arrays */
3996 set_over_alloc_dd(TRUE);
4000 /* PME, if used, is done on all nodes with 1D decomposition */
4002 cr->duty = (DUTY_PP | DUTY_PME);
4003 npme_major = cr->nnodes;
4006 if (inputrec->ePBC == epbcSCREW)
4009 "pbc=%s is only implemented with domain decomposition",
4010 epbc_names[inputrec->ePBC]);
4016 /* After possible communicator splitting in make_dd_communicators.
4017 * we can set up the intra/inter node communication.
4019 gmx_setup_nodecomm(fplog,cr);
4022 wcycle = wallcycle_init(fplog,resetstep,cr);
4025 /* Master synchronizes its value of reset_counters with all nodes
4026 * including PME only nodes */
4027 reset_counters = wcycle_get_reset_counters(wcycle);
4028 gmx_bcast_sim(sizeof(reset_counters),&reset_counters,cr);
4029 wcycle_set_reset_counters(wcycle, reset_counters);
4034 if (cr->duty & DUTY_PP)
4036 /* For domain decomposition we allocate dynamically
4037 * in dd_partition_system.
4039 if (DOMAINDECOMP(cr))
4041 bcast_state_setup(cr,state);
4051 bcast_state(cr,state,TRUE);
4055 /* Dihedral Restraints */
4056 if (gmx_mtop_ftype_count(mtop,F_DIHRES) > 0)
4058 init_dihres(fplog,mtop,inputrec,fcd);
4061 /* Initiate forcerecord */
4063 init_forcerec(fplog,oenv,fr,fcd,inputrec,mtop,cr,box,FALSE,
4064 opt2fn("-table",nfile,fnm),
4065 opt2fn("-tablep",nfile,fnm),
4066 opt2fn("-tableb",nfile,fnm),FALSE,pforce);
4068 /* version for PCA_NOT_READ_NODE (see md.c) */
4069 /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
4070 "nofile","nofile","nofile",FALSE,pforce);
4072 fr->bSepDVDL = ((Flags & MD_SEPPOT) == MD_SEPPOT);
4074 /* Initialize QM-MM */
4077 init_QMMMrec(cr,box,mtop,inputrec,fr);
4080 /* Initialize the mdatoms structure.
4081 * mdatoms is not filled with atom data,
4082 * as this can not be done now with domain decomposition.
4084 mdatoms = init_mdatoms(fplog,mtop,inputrec->efep!=efepNO);
4086 /* Initialize the virtual site communication */
4087 vsite = init_vsite(mtop,cr);
4089 calc_shifts(box,fr->shift_vec);
4091 /* With periodic molecules the charge groups should be whole at start up
4092 * and the virtual sites should not be far from their proper positions.
4094 if (!inputrec->bContinuation && MASTER(cr) &&
4095 !(inputrec->ePBC != epbcNONE && inputrec->bPeriodicMols))
4097 /* Make molecules whole at start of run */
4098 if (fr->ePBC != epbcNONE)
4100 do_pbc_first_mtop(fplog,inputrec->ePBC,box,mtop,state->x);
4104 /* Correct initial vsite positions are required
4105 * for the initial distribution in the domain decomposition
4106 * and for the initial shell prediction.
4108 construct_vsites_mtop(fplog,vsite,mtop,state->x);
4112 /* Initiate PPPM if necessary */
4113 if (fr->eeltype == eelPPPM)
4115 if (mdatoms->nChargePerturbed)
4117 gmx_fatal(FARGS,"Free energy with %s is not implemented",
4118 eel_names[fr->eeltype]);
4120 status = gmx_pppm_init(fplog,cr,oenv,FALSE,TRUE,box,
4121 getenv("GMXGHAT"),inputrec, (Flags & MD_REPRODUCIBLE));
4124 gmx_fatal(FARGS,"Error %d initializing PPPM",status);
4128 if (EEL_PME(fr->eeltype))
4130 ewaldcoeff = fr->ewaldcoeff;
4131 pmedata = &fr->pmedata;
4140 /* This is a PME only node */
4142 /* We don't need the state */
4145 ewaldcoeff = calc_ewaldcoeff(inputrec->rcoulomb, inputrec->ewald_rtol);
4149 /* Initiate PME if necessary,
4150 * either on all nodes or on dedicated PME nodes only. */
4151 if (EEL_PME(inputrec->coulombtype))
4155 nChargePerturbed = mdatoms->nChargePerturbed;
4157 if (cr->npmenodes > 0)
4159 /* The PME only nodes need to know nChargePerturbed */
4160 gmx_bcast_sim(sizeof(nChargePerturbed),&nChargePerturbed,cr);
4162 if (cr->duty & DUTY_PME)
4164 status = gmx_pme_init(pmedata,cr,npme_major,npme_minor,inputrec,
4165 mtop ? mtop->natoms : 0,nChargePerturbed,
4166 (Flags & MD_REPRODUCIBLE));
4169 gmx_fatal(FARGS,"Error %d initializing PME",status);
4175 /* if (integrator[inputrec->eI].func == do_md
4178 integrator[inputrec->eI].func == do_md_openmm
4182 /* Turn on signal handling on all nodes */
4184 * (A user signal from the PME nodes (if any)
4185 * is communicated to the PP nodes.
4187 signal_handler_install();
4190 if (cr->duty & DUTY_PP)
4192 if (inputrec->ePull != epullNO)
4194 /* Initialize pull code */
4195 init_pull(fplog,inputrec,nfile,fnm,mtop,cr,oenv,
4196 EI_DYNAMICS(inputrec->eI) && MASTER(cr),Flags);
4199 constr = init_constraints(fplog,mtop,inputrec,ed,state,cr);
4201 if (DOMAINDECOMP(cr))
4203 dd_init_bondeds(fplog,cr->dd,mtop,vsite,constr,inputrec,
4204 Flags & MD_DDBONDCHECK,fr->cginfo_mb);
4206 set_dd_parameters(fplog,cr->dd,dlb_scale,inputrec,fr,&ddbox);
4208 setup_dd_grid(fplog,cr->dd);
4211 /* Now do whatever the user wants us to do (how flexible...) */
4212 do_md_membed(fplog,cr,nfile,fnm,
4213 oenv,bVerbose,bCompact,
4216 nstepout,inputrec,mtop,
4218 mdatoms,nrnb,wcycle,ed,fr,
4219 repl_ex_nst,repl_ex_seed,
4220 cpt_period,max_hours,
4224 fac, r_ins, pos_ins, ins_at,
4225 xy_step, z_step, it_xy, it_z);
4227 if (inputrec->ePull != epullNO)
4229 finish_pull(fplog,inputrec->pull);
4235 gmx_pmeonly(*pmedata,cr,nrnb,wcycle,ewaldcoeff,FALSE,inputrec);
4238 if (EI_DYNAMICS(inputrec->eI) || EI_TPI(inputrec->eI))
4240 /* Some timing stats */
4243 if (runtime.proc == 0)
4245 runtime.proc = runtime.real;
4254 wallcycle_stop(wcycle,ewcRUN);
4256 /* Finish up, write some stuff
4257 * if rerunMD, don't write last frame again
4259 finish_run(fplog,cr,ftp2fn(efSTO,nfile,fnm),
4260 inputrec,nrnb,wcycle,&runtime,
4261 EI_DYNAMICS(inputrec->eI) && !MULTISIM(cr));
4263 /* Does what it says */
4264 print_date_and_time(fplog,cr->nodeid,"Finished mdrun",&runtime);
4266 /* Close logfile already here if we were appending to it */
4267 if (MASTER(cr) && (Flags & MD_APPENDFILES))
4269 gmx_log_close(fplog);
4277 rc=(int)gmx_get_stop_condition();
4282 int gmx_membed(int argc,char *argv[])
4284 const char *desc[] = {
4285 "[TT]g_membed[tt] embeds a membrane protein into an equilibrated lipid bilayer at the position",
4286 "and orientation specified by the user.\n",
4288 "SHORT MANUAL\n------------\n",
4289 "The user should merge the structure files of the protein and membrane (+solvent), creating a",
4290 "single structure file with the protein overlapping the membrane at the desired position and",
4291 "orientation. Box size should be taken from the membrane structure file. The corresponding topology",
4292 "files should also be merged. Consecutively, create a [TT].tpr[tt] file (input for [TT]g_membed[tt]) from these files,"
4293 "with the following options included in the [TT].mdp[tt] file.\n",
4294 " - [TT]integrator = md[tt][BR]",
4295 " - [TT]energygrp = Protein[tt] (or other group that you want to insert)[BR]",
4296 " - [TT]freezegrps = Protein[tt][BR]",
4297 " - [TT]freezedim = Y Y Y[tt][BR]",
4298 " - [TT]energygrp_excl = Protein Protein[tt][BR]",
4299 "The output is a structure file containing the protein embedded in the membrane. If a topology",
4300 "file is provided, the number of lipid and ",
4301 "solvent molecules will be updated to match the new structure file.\n",
4302 "For a more extensive manual see Wolf et al, J Comp Chem 31 (2010) 2169-2174, Appendix.\n",
4304 "SHORT METHOD DESCRIPTION\n",
4305 "------------------------\n",
4306 "1. The protein is resized around its center of mass by a factor [TT]-xy[tt] in the xy-plane",
4307 "(the membrane plane) and a factor [TT]-z[tt] in the z-direction (if the size of the",
4308 "protein in the z-direction is the same or smaller than the width of the membrane, a",
4309 "[TT]-z[tt] value larger than 1 can prevent that the protein will be enveloped by the lipids).\n",
4310 "2. All lipid and solvent molecules overlapping with the resized protein are removed. All",
4311 "intraprotein interactions are turned off to prevent numerical issues for small values of [TT]-xy[tt]",
4313 "3. One md step is performed.\n",
4314 "4. The resize factor ([TT]-xy[tt] or [TT]-z[tt]) is incremented by a small amount ((1-xy)/nxy or (1-z)/nz) and the",
4315 "protein is resized again around its center of mass. The resize factor for the xy-plane",
4316 "is incremented first. The resize factor for the z-direction is not changed until the [TT]-xy[tt] factor",
4317 "is 1 (thus after [TT]-nxy[tt] iterations).\n",
4318 "5. Repeat step 3 and 4 until the protein reaches its original size ([TT]-nxy[tt] + [TT]-nz[tt] iterations).\n",
4319 "For a more extensive method description see Wolf et al, J Comp Chem, 31 (2010) 2169-2174.\n",
4322 " - Protein can be any molecule you want to insert in the membrane.\n",
4323 " - It is recommended to perform a short equilibration run after the embedding",
4324 "(see Wolf et al, J Comp Chem 31 (2010) 2169-2174), to re-equilibrate the membrane. Clearly",
4325 "protein equilibration might require longer.\n",
4330 { efTPX, "-f", "into_mem", ffREAD },
4331 { efNDX, "-n", "index", ffOPTRD },
4332 { efTOP, "-p", "topol", ffOPTRW },
4333 { efTRN, "-o", NULL, ffWRITE },
4334 { efXTC, "-x", NULL, ffOPTWR },
4335 { efCPT, "-cpi", NULL, ffOPTRD },
4336 { efCPT, "-cpo", NULL, ffOPTWR },
4337 { efSTO, "-c", "membedded", ffWRITE },
4338 { efEDR, "-e", "ener", ffWRITE },
4339 { efLOG, "-g", "md", ffWRITE },
4340 { efEDI, "-ei", "sam", ffOPTRD },
4341 { efTRX, "-rerun", "rerun", ffOPTRD },
4342 { efXVG, "-table", "table", ffOPTRD },
4343 { efXVG, "-tablep", "tablep", ffOPTRD },
4344 { efXVG, "-tableb", "table", ffOPTRD },
4345 { efXVG, "-dhdl", "dhdl", ffOPTWR },
4346 { efXVG, "-field", "field", ffOPTWR },
4347 { efXVG, "-table", "table", ffOPTRD },
4348 { efXVG, "-tablep", "tablep", ffOPTRD },
4349 { efXVG, "-tableb", "table", ffOPTRD },
4350 { efTRX, "-rerun", "rerun", ffOPTRD },
4351 { efXVG, "-tpi", "tpi", ffOPTWR },
4352 { efXVG, "-tpid", "tpidist", ffOPTWR },
4353 { efEDI, "-ei", "sam", ffOPTRD },
4354 { efEDO, "-eo", "sam", ffOPTWR },
4355 { efGCT, "-j", "wham", ffOPTRD },
4356 { efGCT, "-jo", "bam", ffOPTWR },
4357 { efXVG, "-ffout", "gct", ffOPTWR },
4358 { efXVG, "-devout", "deviatie", ffOPTWR },
4359 { efXVG, "-runav", "runaver", ffOPTWR },
4360 { efXVG, "-px", "pullx", ffOPTWR },
4361 { efXVG, "-pf", "pullf", ffOPTWR },
4362 { efMTX, "-mtx", "nm", ffOPTWR },
4363 { efNDX, "-dn", "dipole", ffOPTWR },
4364 { efRND, "-multidir",NULL, ffOPTRDMULT}
4366 #define NFILE asize(fnm)
4368 /* Command line options ! */
4369 gmx_bool bCart = FALSE;
4370 gmx_bool bPPPME = FALSE;
4371 gmx_bool bPartDec = FALSE;
4372 gmx_bool bDDBondCheck = TRUE;
4373 gmx_bool bDDBondComm = TRUE;
4374 gmx_bool bVerbose = FALSE;
4375 gmx_bool bCompact = TRUE;
4376 gmx_bool bSepPot = FALSE;
4377 gmx_bool bRerunVSite = FALSE;
4378 gmx_bool bIonize = FALSE;
4379 gmx_bool bConfout = TRUE;
4380 gmx_bool bReproducible = FALSE;
4384 int nstglobalcomm=-1;
4386 int repl_ex_seed=-1;
4388 int nthreads=0; /* set to determine # of threads automatically */
4391 rvec realddxyz={0,0,0};
4392 const char *ddno_opt[ddnoNR+1] =
4393 { NULL, "interleave", "pp_pme", "cartesian", NULL };
4394 const char *dddlb_opt[] =
4395 { NULL, "auto", "no", "yes", NULL };
4396 real rdd=0.0,rconstr=0.0,dlb_scale=0.8,pforce=-1;
4397 char *ddcsx=NULL,*ddcsy=NULL,*ddcsz=NULL;
4398 real cpt_period=15.0,max_hours=-1;
4399 gmx_bool bAppendFiles=TRUE,bAddPart=TRUE;
4400 gmx_bool bResetCountersHalfWay=FALSE;
4401 output_env_t oenv=NULL;
4402 const char *deviceOptions = "";
4410 real probe_rad = 0.22;
4414 gmx_bool bALLOW_ASYMMETRY=FALSE;
4417 /* arguments relevant to OPENMM only*/
4419 gmx_input("g_membed not functional in openmm");
4423 { "-xyinit", FALSE, etREAL, {&xy_fac}, "Resize factor for the protein in the xy dimension before starting embedding" },
4424 { "-xyend", FALSE, etREAL, {&xy_max}, "Final resize factor in the xy dimension" },
4425 { "-zinit", FALSE, etREAL, {&z_fac}, "Resize factor for the protein in the z dimension before starting embedding" },
4426 { "-zend", FALSE, etREAL, {&z_max}, "Final resize faction in the z dimension" },
4427 { "-nxy", FALSE, etINT, {&it_xy}, "Number of iteration for the xy dimension" },
4428 { "-nz", FALSE, etINT, {&it_z}, "Number of iterations for the z dimension" },
4429 { "-rad", FALSE, etREAL, {&probe_rad}, "Probe radius to check for overlap between the group to embed and the membrane"},
4430 { "-pieces", FALSE, etINT, {&pieces}, "Perform piecewise resize. Select parts of the group to insert and resize these with respect to their own geometrical center." },
4431 { "-asymmetry",FALSE, etBOOL,{&bALLOW_ASYMMETRY}, "Allow asymmetric insertion, i.e. the number of lipids removed from the upper and lower leaflet will not be checked." },
4432 { "-ndiff" , FALSE, etINT, {&low_up_rm}, "Number of lipids that will additionally be removed from the lower (negative number) or upper (positive number) membrane leaflet." },
4433 { "-maxwarn", FALSE, etINT, {&maxwarn}, "Maximum number of warning allowed" },
4434 { "-pd", FALSE, etBOOL,{&bPartDec},
4435 "HIDDENUse particle decompostion" },
4436 { "-dd", FALSE, etRVEC,{&realddxyz},
4437 "HIDDENDomain decomposition grid, 0 is optimize" },
4438 { "-nt", FALSE, etINT, {&nthreads},
4439 "HIDDENNumber of threads to start (0 is guess)" },
4440 { "-npme", FALSE, etINT, {&npme},
4441 "HIDDENNumber of separate nodes to be used for PME, -1 is guess" },
4442 { "-ddorder", FALSE, etENUM, {ddno_opt},
4443 "HIDDENDD node order" },
4444 { "-ddcheck", FALSE, etBOOL, {&bDDBondCheck},
4445 "HIDDENCheck for all bonded interactions with DD" },
4446 { "-ddbondcomm", FALSE, etBOOL, {&bDDBondComm},
4447 "HIDDENUse special bonded atom communication when [TT]-rdd[tt] > cut-off" },
4448 { "-rdd", FALSE, etREAL, {&rdd},
4449 "HIDDENThe maximum distance for bonded interactions with DD (nm), 0 is determine from initial coordinates" },
4450 { "-rcon", FALSE, etREAL, {&rconstr},
4451 "HIDDENMaximum distance for P-LINCS (nm), 0 is estimate" },
4452 { "-dlb", FALSE, etENUM, {dddlb_opt},
4453 "HIDDENDynamic load balancing (with DD)" },
4454 { "-dds", FALSE, etREAL, {&dlb_scale},
4455 "HIDDENMinimum allowed dlb scaling of the DD cell size" },
4456 { "-ddcsx", FALSE, etSTR, {&ddcsx},
4457 "HIDDENThe DD cell sizes in x" },
4458 { "-ddcsy", FALSE, etSTR, {&ddcsy},
4459 "HIDDENThe DD cell sizes in y" },
4460 { "-ddcsz", FALSE, etSTR, {&ddcsz},
4461 "HIDDENThe DD cell sizes in z" },
4462 { "-gcom", FALSE, etINT,{&nstglobalcomm},
4463 "HIDDENGlobal communication frequency" },
4464 { "-compact", FALSE, etBOOL,{&bCompact},
4465 "Write a compact log file" },
4466 { "-seppot", FALSE, etBOOL, {&bSepPot},
4467 "HIDDENWrite separate V and dVdl terms for each interaction type and node to the log file(s)" },
4468 { "-pforce", FALSE, etREAL, {&pforce},
4469 "HIDDENPrint all forces larger than this (kJ/mol nm)" },
4470 { "-reprod", FALSE, etBOOL,{&bReproducible},
4471 "HIDDENTry to avoid optimizations that affect binary reproducibility" },
4472 { "-multi", FALSE, etINT,{&nmultisim},
4473 "HIDDENDo multiple simulations in parallel" },
4474 { "-replex", FALSE, etINT, {&repl_ex_nst},
4475 "HIDDENAttempt replica exchange every # steps" },
4476 { "-reseed", FALSE, etINT, {&repl_ex_seed},
4477 "HIDDENSeed for replica exchange, -1 is generate a seed" },
4478 { "-rerunvsite", FALSE, etBOOL, {&bRerunVSite},
4479 "HIDDENRecalculate virtual site coordinates with [TT]-rerun[tt]" },
4480 { "-ionize", FALSE, etBOOL,{&bIonize},
4481 "HIDDENDo a simulation including the effect of an X-Ray bombardment on your system" },
4482 { "-confout", TRUE, etBOOL, {&bConfout},
4483 "HIDDENWrite the last configuration with [TT]-c[tt] and force checkpointing at the last step" },
4484 { "-stepout", FALSE, etINT, {&nstepout},
4485 "HIDDENFrequency of writing the remaining runtime" },
4486 { "-resetstep", FALSE, etINT, {&resetstep},
4487 "HIDDENReset cycle counters after these many time steps" },
4488 { "-resethway", FALSE, etBOOL, {&bResetCountersHalfWay},
4489 "HIDDENReset the cycle counters after half the number of steps or halfway [TT]-maxh[tt]" },
4490 { "-v", FALSE, etBOOL,{&bVerbose},
4491 "Be loud and noisy" },
4492 { "-maxh", FALSE, etREAL, {&max_hours},
4493 "HIDDENTerminate after 0.99 times this time (hours)" },
4494 { "-cpt", FALSE, etREAL, {&cpt_period},
4495 "HIDDENCheckpoint interval (minutes)" },
4496 { "-append", FALSE, etBOOL, {&bAppendFiles},
4497 "HIDDENAppend to previous output files when continuing from checkpoint" },
4498 { "-addpart", FALSE, etBOOL, {&bAddPart},
4499 "HIDDENAdd the simulation part number to all output files when continuing from checkpoint" },
4502 unsigned long Flags, PCA_Flags;
4505 gmx_bool HaveCheckpoint;
4506 FILE *fplog,*fptest;
4507 int sim_part,sim_part_fn;
4508 const char *part_suffix=".part";
4509 char suffix[STRLEN];
4511 char **multidir=NULL;
4513 cr = init_par(&argc,&argv);
4515 PCA_Flags = (PCA_KEEP_ARGS | PCA_NOEXIT_ON_ARGS | PCA_CAN_SET_DEFFNM
4516 | (MASTER(cr) ? 0 : PCA_QUIET));
4519 /* Comment this in to do fexist calls only on master
4520 * works not with rerun or tables at the moment
4521 * also comment out the version of init_forcerec in md.c
4522 * with NULL instead of opt2fn
4527 PCA_Flags |= PCA_NOT_READ_NODE;
4531 parse_common_args(&argc,argv,PCA_Flags, NFILE,fnm,asize(pa),pa,
4532 asize(desc),desc,0,NULL, &oenv);
4534 /* we set these early because they might be used in init_multisystem()
4535 Note that there is the potential for npme>nnodes until the number of
4536 threads is set later on, if there's thread parallelization. That shouldn't
4537 lead to problems. */
4538 dd_node_order = nenum(ddno_opt);
4539 cr->npmenodes = npme;
4542 /* now determine the number of threads automatically. The threads are
4543 only started at mdrunner_threads, though. */
4546 nthreads=tMPI_Get_recommended_nthreads();
4552 /* now check the -multi and -multidir option */
4553 if (opt2bSet("-multidir", NFILE, fnm))
4558 gmx_fatal(FARGS, "mdrun -multi and -multidir options are mutually exclusive.");
4560 nmultisim = opt2fns(&multidir, "-multidir", NFILE, fnm);
4564 if (repl_ex_nst != 0 && nmultisim < 2)
4565 gmx_fatal(FARGS,"Need at least two replicas for replica exchange (option -multi)");
4567 if (nmultisim > 1) {
4569 gmx_bool bParFn = (multidir == NULL);
4570 init_multisystem(cr,nmultisim,multidir,NFILE,fnm,TRUE);
4572 gmx_fatal(FARGS,"mdrun -multi is not supported with the thread library.Please compile GROMACS with MPI support");
4576 /* Check if there is ANY checkpoint file available */
4578 sim_part_fn = sim_part;
4579 if (opt2bSet("-cpi",NFILE,fnm))
4582 read_checkpoint_simulation_part(opt2fn_master("-cpi", NFILE,fnm,cr),
4583 &sim_part_fn,NULL,cr,
4584 bAppendFiles,NFILE,fnm,
4585 part_suffix,&bAddPart);
4586 if (sim_part_fn==0 && MASTER(cr))
4588 fprintf(stdout,"No previous checkpoint file present, assuming this is a new run.\n");
4592 sim_part = sim_part_fn + 1;
4597 bAppendFiles = FALSE;
4602 sim_part_fn = sim_part;
4605 if (bAddPart && sim_part_fn > 1)
4607 /* This is a continuation run, rename trajectory output files
4608 (except checkpoint files) */
4609 /* create new part name first (zero-filled) */
4610 sprintf(suffix,"%s%04d",part_suffix,sim_part_fn);
4612 add_suffix_to_output_names(fnm,NFILE,suffix);
4613 fprintf(stdout,"Checkpoint file is from part %d, new output files will be suffixed '%s'.\n",sim_part-1,suffix);
4616 Flags = opt2bSet("-rerun",NFILE,fnm) ? MD_RERUN : 0;
4617 Flags = Flags | (bSepPot ? MD_SEPPOT : 0);
4618 Flags = Flags | (bIonize ? MD_IONIZE : 0);
4619 Flags = Flags | (bPartDec ? MD_PARTDEC : 0);
4620 Flags = Flags | (bDDBondCheck ? MD_DDBONDCHECK : 0);
4621 Flags = Flags | (bDDBondComm ? MD_DDBONDCOMM : 0);
4622 Flags = Flags | (bConfout ? MD_CONFOUT : 0);
4623 Flags = Flags | (bRerunVSite ? MD_RERUN_VSITE : 0);
4624 Flags = Flags | (bReproducible ? MD_REPRODUCIBLE : 0);
4625 Flags = Flags | (bAppendFiles ? MD_APPENDFILES : 0);
4626 Flags = Flags | (sim_part>1 ? MD_STARTFROMCPT : 0);
4627 Flags = Flags | (bResetCountersHalfWay ? MD_RESETCOUNTERSHALFWAY : 0);
4630 /* We postpone opening the log file if we are appending, so we can
4631 first truncate the old log file and append to the correct position
4633 if ((MASTER(cr) || bSepPot) && !bAppendFiles)
4635 gmx_log_open(ftp2fn(efLOG,NFILE,fnm),cr,!bSepPot,Flags,&fplog);
4636 CopyRight(fplog,argv[0]);
4637 please_cite(fplog,"Hess2008b");
4638 please_cite(fplog,"Spoel2005a");
4639 please_cite(fplog,"Lindahl2001a");
4640 please_cite(fplog,"Berendsen95a");
4647 ddxyz[XX] = (int)(realddxyz[XX] + 0.5);
4648 ddxyz[YY] = (int)(realddxyz[YY] + 0.5);
4649 ddxyz[ZZ] = (int)(realddxyz[ZZ] + 0.5);
4651 /* even if nthreads = 1, we still call this one */
4653 rc = mdrunner_membed(fplog, cr, NFILE, fnm, oenv, bVerbose, bCompact,
4655 ddxyz, dd_node_order, rdd, rconstr, dddlb_opt[0], dlb_scale,
4656 ddcsx, ddcsy, ddcsz, nstepout, resetstep, nmultisim, repl_ex_nst,
4657 repl_ex_seed, pforce, cpt_period, max_hours, deviceOptions, Flags,
4658 xy_fac,xy_max,z_fac,z_max,
4659 it_xy,it_z,probe_rad,low_up_rm,
4660 pieces,bALLOW_ASYMMETRY,maxwarn);
4662 if (gmx_parallel_env_initialized())
4665 if (MULTIMASTER(cr)) {
4669 /* Log file has to be closed in mdrunner if we are appending to it
4670 (fplog not set here) */
4671 fprintf(stderr,"Please cite:\nWolf et al, J Comp Chem 31 (2010) 2169-2174.\n");
4673 if (MASTER(cr) && !bAppendFiles)
4675 gmx_log_close(fplog);