2 * $Id: mdrun.c,v 1.139.2.9 2009/05/04 16:13:29 hess Exp $
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
52 #include "checkpoint.h"
77 #include "mpelogging.h"
85 #include "checkpoint.h"
86 #include "mtop_util.h"
89 #include "sighandler.h"
102 /* We use the same defines as in mvdata.c here */
103 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d),(cr))
104 #define nblock_bc(cr,nr,d) gmx_bcast((nr)*sizeof((d)[0]), (d),(cr))
105 #define snew_bc(cr,d,nr) { if (!MASTER(cr)) snew((d),(nr)); }
107 /* The following two variables and the signal_handler function
108 * is used from pme.c as well
137 int natoms
; /*nr of atoms per lipid*/
138 int mol1
; /*id of the first lipid molecule*/
159 int search_string(char *s
,int ng
,char ***gn
)
163 for(i
=0; (i
<ng
); i
++)
164 if (gmx_strcasecmp(s
,*gn
[i
]) == 0)
167 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
);
172 int get_mol_id(int at
,int nmblock
,gmx_molblock_t
*mblock
, int *type
, int *block
)
177 for(i
=0;i
<nmblock
;i
++)
179 if(at
<(mblock
[i
].nmol
*mblock
[i
].natoms_mol
))
181 mol_id
+=at
/mblock
[i
].natoms_mol
;
182 *type
= mblock
[i
].type
;
186 at
-= mblock
[i
].nmol
*mblock
[i
].natoms_mol
;
187 mol_id
+=mblock
[i
].nmol
;
191 gmx_fatal(FARGS
,"Something is wrong in mol ids, at %d, mol_id %d",at
,mol_id
);
196 int get_block(int mol_id
,int nmblock
,gmx_molblock_t
*mblock
)
201 for(i
=0;i
<nmblock
;i
++)
203 nmol
+=mblock
[i
].nmol
;
208 gmx_fatal(FARGS
,"mol_id %d larger than total number of molecules %d.\n",mol_id
,nmol
);
213 int get_tpr_version(const char *infile
)
220 fio
= open_tpx(infile
,"r");
221 gmx_fio_checktype(fio
);
223 precision
= sizeof(real
);
225 gmx_fio_do_string(fio
,buf
);
226 if (strncmp(buf
,"VERSION",7))
227 gmx_fatal(FARGS
,"Can not read file %s,\n"
228 " this file is from a Gromacs version which is older than 2.0\n"
229 " Make a new one with grompp or use a gro or pdb file, if possible",
230 gmx_fio_getname(fio
));
231 gmx_fio_do_int(fio
,precision
);
232 bDouble
= (precision
== sizeof(double));
233 if ((precision
!= sizeof(float)) && !bDouble
)
234 gmx_fatal(FARGS
,"Unknown precision in file %s: real is %d bytes "
235 "instead of %d or %d",
236 gmx_fio_getname(fio
),precision
,sizeof(float),sizeof(double));
237 gmx_fio_setprecision(fio
,bDouble
);
238 fprintf(stderr
,"Reading file %s, %s (%s precision)\n",
239 gmx_fio_getname(fio
),buf
,bDouble
? "double" : "single");
241 gmx_fio_do_int(fio
,fver
);
248 void set_inbox(int natom
, rvec
*x
)
253 tmp
[XX
]=tmp
[YY
]=tmp
[ZZ
]=0.0;
256 if(x
[i
][XX
]<tmp
[XX
]) tmp
[XX
]=x
[i
][XX
];
257 if(x
[i
][YY
]<tmp
[YY
]) tmp
[YY
]=x
[i
][YY
];
258 if(x
[i
][ZZ
]<tmp
[ZZ
]) tmp
[ZZ
]=x
[i
][ZZ
];
265 int get_mtype_list(t_block
*at
, gmx_mtop_t
*mtop
, t_block
*tlist
)
272 snew(tlist
->index
,at
->nr
);
273 for (i
=0;i
<at
->nr
;i
++)
276 mol_id
= get_mol_id(at
->index
[i
],mtop
->nmolblock
,mtop
->molblock
,&type
,&block
);
279 if(tlist
->index
[j
]==type
)
284 tlist
->index
[nr
]=type
;
289 srenew(tlist
->index
,nr
);
293 void check_types(t_block
*ins_at
,t_block
*rest_at
,gmx_mtop_t
*mtop
)
295 t_block
*ins_mtype
,*rest_mtype
;
300 ins_mtype
->nr
= get_mtype_list(ins_at
, mtop
, ins_mtype
);
301 rest_mtype
->nr
= get_mtype_list(rest_at
, mtop
, rest_mtype
);
303 for(i
=0;i
<ins_mtype
->nr
;i
++)
305 for(j
=0;j
<rest_mtype
->nr
;j
++)
307 if(ins_mtype
->index
[i
]==rest_mtype
->index
[j
])
308 gmx_fatal(FARGS
,"Moleculetype %s is found both in the group to insert and the rest of the system.\n"
309 "Because we need to exclude all interactions between the atoms in the group to\n"
310 "insert, the same moleculetype can not be used in both groups. Change the\n"
311 "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n"
312 "an appropriate *.itp file",*(mtop
->moltype
[rest_mtype
->index
[j
]].name
),
313 *(mtop
->moltype
[rest_mtype
->index
[j
]].name
));
317 sfree(ins_mtype
->index
);
318 sfree(rest_mtype
->index
);
323 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
)
326 real x
,xmin
,xmax
,y
,ymin
,ymax
,z
,zmin
,zmax
;
328 snew(rest_at
->index
,state
->natoms
);
330 xmin
=xmax
=state
->x
[ins_at
->index
[0]][XX
];
331 ymin
=ymax
=state
->x
[ins_at
->index
[0]][YY
];
332 zmin
=zmax
=state
->x
[ins_at
->index
[0]][ZZ
];
334 for(i
=0;i
<state
->natoms
;i
++)
336 gid
= groups
->grpnr
[egcFREEZE
][i
];
337 if(groups
->grps
[egcFREEZE
].nm_ind
[gid
]==ins_grp_id
)
355 srenew(rest_at
->index
,c
);
359 pos_ins
->xmin
[XX
]=xmin
-((xmax
-xmin
)*xy_max
-(xmax
-xmin
))/2;
360 pos_ins
->xmin
[YY
]=ymin
-((ymax
-ymin
)*xy_max
-(ymax
-ymin
))/2;
362 pos_ins
->xmax
[XX
]=xmax
+((xmax
-xmin
)*xy_max
-(xmax
-xmin
))/2;
363 pos_ins
->xmax
[YY
]=ymax
+((ymax
-ymin
)*xy_max
-(ymax
-ymin
))/2;
365 pos_ins
->xmin
[XX
]=xmin
;
366 pos_ins
->xmin
[YY
]=ymin
;
368 pos_ins
->xmax
[XX
]=xmax
;
369 pos_ins
->xmax
[YY
]=ymax
;
372 /* 6.0 is estimated thickness of bilayer */
373 if( (zmax
-zmin
) < 6.0 )
375 pos_ins
->xmin
[ZZ
]=zmin
+(zmax
-zmin
)/2.0-3.0;
376 pos_ins
->xmax
[ZZ
]=zmin
+(zmax
-zmin
)/2.0+3.0;
378 pos_ins
->xmin
[ZZ
]=zmin
;
379 pos_ins
->xmax
[ZZ
]=zmax
;
385 real
est_prot_area(pos_ins_t
*pos_ins
,rvec
*r
,t_block
*ins_at
, mem_t
*mem_p
)
387 real x
,y
,dx
=0.15,dy
=0.15,area
=0.0;
391 for(x
=pos_ins
->xmin
[XX
];x
<pos_ins
->xmax
[XX
];x
+=dx
)
393 for(y
=pos_ins
->xmin
[YY
];y
<pos_ins
->xmax
[YY
];y
+=dy
)
400 if ( (r
[at
][XX
]>=x
) && (r
[at
][XX
]<x
+dx
) &&
401 (r
[at
][YY
]>=y
) && (r
[at
][YY
]<y
+dy
) &&
402 (r
[at
][ZZ
]>mem_p
->zmin
+1.0) && (r
[at
][ZZ
]<mem_p
->zmax
-1.0) )
405 } while ( (c
<ins_at
->nr
) && (add
<0.5) );
414 void init_lip(matrix box
, gmx_mtop_t
*mtop
, lip_t
*lip
)
420 mem_area
= box
[XX
][XX
]*box
[YY
][YY
]-box
[XX
][YY
]*box
[YY
][XX
];
421 for(i
=0;i
<mtop
->nmolblock
;i
++)
423 if(mtop
->molblock
[i
].type
== lip
->id
)
425 lip
->nr
=mtop
->molblock
[i
].nmol
;
426 lip
->natoms
=mtop
->molblock
[i
].natoms_mol
;
429 lip
->area
=2.0*mem_area
/(double)lip
->nr
;
431 for (i
=0;i
<lip
->id
;i
++)
432 mol1
+=mtop
->molblock
[i
].nmol
;
436 int init_mem_at(mem_t
*mem_p
, gmx_mtop_t
*mtop
, rvec
*r
, matrix box
, pos_ins_t
*pos_ins
)
438 int i
,j
,at
,mol
,nmol
,nmolbox
,count
;
440 real z
,zmin
,zmax
,mem_area
;
446 mem_a
=&(mem_p
->mem_at
);
447 snew(mol_id
,mem_a
->nr
);
448 /* snew(index,mem_a->nr); */
449 zmin
=pos_ins
->xmax
[ZZ
];
450 zmax
=pos_ins
->xmin
[ZZ
];
451 for(i
=0;i
<mem_a
->nr
;i
++)
454 if( (r
[at
][XX
]>pos_ins
->xmin
[XX
]) && (r
[at
][XX
]<pos_ins
->xmax
[XX
]) &&
455 (r
[at
][YY
]>pos_ins
->xmin
[YY
]) && (r
[at
][YY
]<pos_ins
->xmax
[YY
]) &&
456 (r
[at
][ZZ
]>pos_ins
->xmin
[ZZ
]) && (r
[at
][ZZ
]<pos_ins
->xmax
[ZZ
]) )
458 mol
= get_mol_id(at
,mtop
->nmolblock
,mtop
->molblock
,&type
,&block
);
475 /* index[count]=at;*/
482 mem_p
->mol_id
=mol_id
;
483 /* srenew(index,count);*/
484 /* mem_p->mem_at.nr=count;*/
485 /* sfree(mem_p->mem_at.index);*/
486 /* mem_p->mem_at.index=index;*/
488 if((zmax
-zmin
)>(box
[ZZ
][ZZ
]-0.5))
489 gmx_fatal(FARGS
,"Something is wrong with your membrane. Max and min z values are %f and %f.\n"
490 "Maybe your membrane is not centered in the box, but located at the box edge in the z-direction,\n"
491 "so that one membrane is distributed over two periodic box images. Another possibility is that\n"
492 "your water layer is not thick enough.\n",zmax
,zmin
);
495 mem_p
->zmed
=(zmax
-zmin
)/2+zmin
;
497 /*number of membrane molecules in protein box*/
498 nmolbox
= count
/mtop
->molblock
[block
].natoms_mol
;
499 /*mem_area = box[XX][XX]*box[YY][YY]-box[XX][YY]*box[YY][XX];
500 mem_p->lip_area = 2.0*mem_area/(double)mem_p->nmol;*/
501 mem_area
= (pos_ins
->xmax
[XX
]-pos_ins
->xmin
[XX
])*(pos_ins
->xmax
[YY
]-pos_ins
->xmin
[YY
]);
502 mem_p
->lip_area
= 2.0*mem_area
/(double)nmolbox
;
504 return mem_p
->mem_at
.nr
;
507 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
)
509 int i
,j
,at
,c
,outsidesum
,gctr
=0;
513 for (i
=0;i
<pos_ins
->pieces
;i
++)
514 idxsum
+=pos_ins
->nidx
[i
];
515 if (idxsum
!=ins_at
->nr
)
516 gmx_fatal(FARGS
,"Piecewise sum of inserted atoms not same as size of group selected to insert.");
518 snew(pos_ins
->geom_cent
,pos_ins
->pieces
);
519 for (i
=0;i
<pos_ins
->pieces
;i
++)
524 pos_ins
->geom_cent
[i
][j
]=0;
527 pos_ins
->geom_cent
[i
][j
]=0;
528 for (j
=0;j
<pos_ins
->nidx
[i
];j
++)
530 at
=pos_ins
->subindex
[i
][j
];
531 copy_rvec(r
[at
],r_ins
[gctr
]);
532 if( (r_ins
[gctr
][ZZ
]<mem_p
->zmax
) && (r_ins
[gctr
][ZZ
]>mem_p
->zmin
) )
534 rvec_inc(pos_ins
->geom_cent
[i
],r_ins
[gctr
]);
542 svmul(1/(double)c
,pos_ins
->geom_cent
[i
],pos_ins
->geom_cent
[i
]);
543 if (!bALLOW_ASYMMETRY
)
544 pos_ins
->geom_cent
[i
][ZZ
]=mem_p
->zmed
;
546 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
]);
548 fprintf(stderr
,"\n");
551 void resize(t_block
*ins_at
, rvec
*r_ins
, rvec
*r
, pos_ins_t
*pos_ins
,rvec fac
)
554 for (k
=0;k
<pos_ins
->pieces
;k
++)
555 for(i
=0;i
<pos_ins
->nidx
[k
];i
++)
557 at
=pos_ins
->subindex
[k
][i
];
559 r
[at
][j
]=pos_ins
->geom_cent
[k
][j
]+fac
[j
]*(r_ins
[c
][j
]-pos_ins
->geom_cent
[k
][j
]);
564 int gen_rm_list(rm_t
*rm_p
,t_block
*ins_at
,t_block
*rest_at
,t_pbc
*pbc
, gmx_mtop_t
*mtop
,
565 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
)
567 int i
,j
,k
,l
,at
,at2
,mol_id
;
569 int nrm
,nupper
,nlower
;
570 real r_min_rad
,z_lip
,min_norm
;
576 r_min_rad
=probe_rad
*probe_rad
;
577 snew(rm_p
->mol
,mtop
->mols
.nr
);
578 snew(rm_p
->block
,mtop
->mols
.nr
);
581 for(i
=0;i
<ins_at
->nr
;i
++)
584 for(j
=0;j
<rest_at
->nr
;j
++)
586 at2
=rest_at
->index
[j
];
587 pbc_dx(pbc
,r
[at
],r
[at2
],dr
);
589 if(norm2(dr
)<r_min_rad
)
591 mol_id
= get_mol_id(at2
,mtop
->nmolblock
,mtop
->molblock
,&type
,&block
);
594 if(rm_p
->mol
[l
]==mol_id
)
598 /*fprintf(stderr,"%d wordt toegevoegd\n",mol_id);*/
599 rm_p
->mol
[nrm
]=mol_id
;
600 rm_p
->block
[nrm
]=block
;
603 for(l
=0;l
<mem_p
->nmol
;l
++)
605 if(mol_id
==mem_p
->mol_id
[l
])
607 for(k
=mtop
->mols
.index
[mol_id
];k
<mtop
->mols
.index
[mol_id
+1];k
++)
609 z_lip
/=mtop
->molblock
[block
].natoms_mol
;
610 if(z_lip
<mem_p
->zmed
)
621 /*make sure equal number of lipids from upper and lower layer are removed */
622 if( (nupper
!=nlower
) && (!bALLOW_ASYMMETRY
) )
624 snew(dist
,mem_p
->nmol
);
625 snew(order
,mem_p
->nmol
);
626 for(i
=0;i
<mem_p
->nmol
;i
++)
628 at
= mtop
->mols
.index
[mem_p
->mol_id
[i
]];
629 pbc_dx(pbc
,r
[at
],pos_ins
->geom_cent
[0],dr
);
630 if (pos_ins
->pieces
>1)
634 for (k
=1;k
<pos_ins
->pieces
;k
++)
636 pbc_dx(pbc
,r
[at
],pos_ins
->geom_cent
[k
],dr_tmp
);
637 if (norm2(dr_tmp
) < min_norm
)
639 min_norm
=norm2(dr_tmp
);
640 copy_rvec(dr_tmp
,dr
);
644 dist
[i
]=dr
[XX
]*dr
[XX
]+dr
[YY
]*dr
[YY
];
646 while (j
>=0 && dist
[i
]<dist
[order
[j
]])
655 while(nupper
!=nlower
)
657 mol_id
=mem_p
->mol_id
[order
[i
]];
658 block
=get_block(mol_id
,mtop
->nmolblock
,mtop
->molblock
);
662 if(rm_p
->mol
[l
]==mol_id
)
667 for(k
=mtop
->mols
.index
[mol_id
];k
<mtop
->mols
.index
[mol_id
+1];k
++)
669 z_lip
/=mtop
->molblock
[block
].natoms_mol
;
670 if(nupper
>nlower
&& z_lip
<mem_p
->zmed
)
672 rm_p
->mol
[nrm
]=mol_id
;
673 rm_p
->block
[nrm
]=block
;
677 else if (nupper
<nlower
&& z_lip
>mem_p
->zmed
)
679 rm_p
->mol
[nrm
]=mol_id
;
680 rm_p
->block
[nrm
]=block
;
688 gmx_fatal(FARGS
,"Trying to remove more lipid molecules than there are in the membrane");
695 srenew(rm_p
->mol
,nrm
);
696 srenew(rm_p
->block
,nrm
);
698 return nupper
+nlower
;
701 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
)
703 int i
,j
,k
,n
,rm
,mol_id
,at
,block
;
705 atom_id
*list
,*new_mols
;
706 unsigned char *new_egrp
[egcNR
];
709 snew(list
,state
->natoms
);
711 for(i
=0;i
<rm_p
->nr
;i
++)
714 at
=mtop
->mols
.index
[mol_id
];
715 block
=rm_p
->block
[i
];
716 mtop
->molblock
[block
].nmol
--;
717 for(j
=0;j
<mtop
->molblock
[block
].natoms_mol
;j
++)
723 mtop
->mols
.index
[mol_id
]=-1;
726 mtop
->mols
.nr
-=rm_p
->nr
;
727 mtop
->mols
.nalloc_index
-=rm_p
->nr
;
728 snew(new_mols
,mtop
->mols
.nr
);
729 for(i
=0;i
<mtop
->mols
.nr
+rm_p
->nr
;i
++)
732 if(mtop
->mols
.index
[i
]!=-1)
734 new_mols
[j
]=mtop
->mols
.index
[i
];
738 sfree(mtop
->mols
.index
);
739 mtop
->mols
.index
=new_mols
;
744 state
->nalloc
=state
->natoms
;
745 snew(x_tmp
,state
->nalloc
);
746 snew(v_tmp
,state
->nalloc
);
750 if(groups
->grpnr
[i
]!=NULL
)
752 groups
->ngrpnr
[i
]=state
->natoms
;
753 snew(new_egrp
[i
],state
->natoms
);
758 for (i
=0;i
<state
->natoms
+n
;i
++)
774 if(groups
->grpnr
[j
]!=NULL
)
776 new_egrp
[j
][i
-rm
]=groups
->grpnr
[j
][i
];
779 copy_rvec(state
->x
[i
],x_tmp
[i
-rm
]);
780 copy_rvec(state
->v
[i
],v_tmp
[i
-rm
]);
781 for(j
=0;j
<ins_at
->nr
;j
++)
783 if (i
==ins_at
->index
[j
])
784 ins_at
->index
[j
]=i
-rm
;
786 for(j
=0;j
<pos_ins
->pieces
;j
++)
788 for(k
=0;k
<pos_ins
->nidx
[j
];k
++)
790 if (i
==pos_ins
->subindex
[j
][k
])
791 pos_ins
->subindex
[j
][k
]=i
-rm
;
803 if(groups
->grpnr
[i
]!=NULL
)
805 sfree(groups
->grpnr
[i
]);
806 groups
->grpnr
[i
]=new_egrp
[i
];
811 int rm_bonded(t_block
*ins_at
, gmx_mtop_t
*mtop
)
814 int type
,natom
,nmol
,at
,atom1
=0,rm_at
=0;
816 /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/
817 /*this routine does not live as dangerously as it seems. There is namely a check in mdrunner_membed to make
818 *sure that g_membed exits with a warning when there are molecules of the same type not in the
819 *ins_at index group. MGWolf 050710 */
822 snew(bRM
,mtop
->nmoltype
);
823 for (i
=0;i
<mtop
->nmoltype
;i
++)
828 for (i
=0;i
<mtop
->nmolblock
;i
++)
830 /*loop over molecule blocks*/
831 type
=mtop
->molblock
[i
].type
;
832 natom
=mtop
->molblock
[i
].natoms_mol
;
833 nmol
=mtop
->molblock
[i
].nmol
;
835 for(j
=0;j
<natom
*nmol
&& bRM
[type
]==TRUE
;j
++)
837 /*loop over atoms in the block*/
838 at
=j
+atom1
; /*atom index = block index + offset*/
841 for (m
=0;(m
<ins_at
->nr
) && (bINS
==FALSE
);m
++)
843 /*loop over atoms in insertion index group to determine if we're inserting one*/
844 if(at
==ins_at
->index
[m
])
851 atom1
+=natom
*nmol
; /*update offset*/
854 rm_at
+=natom
*nmol
; /*increment bonded removal counter by # atoms in block*/
858 for(i
=0;i
<mtop
->nmoltype
;i
++)
864 mtop
->moltype
[i
].ilist
[j
].nr
=0;
866 for(j
=F_POSRES
;j
<=F_VSITEN
;j
++)
868 mtop
->moltype
[i
].ilist
[j
].nr
=0;
877 void top_update(const char *topfile
, char *ins
, rm_t
*rm_p
, gmx_mtop_t
*mtop
)
879 #define TEMP_FILENM "temp.top"
882 char buf
[STRLEN
],buf2
[STRLEN
],*temp
;
883 int i
,*nmol_rm
,nmol
,line
;
885 fpin
= ffopen(topfile
,"r");
886 fpout
= ffopen(TEMP_FILENM
,"w");
888 snew(nmol_rm
,mtop
->nmoltype
);
889 for(i
=0;i
<rm_p
->nr
;i
++)
890 nmol_rm
[rm_p
->block
[i
]]++;
893 while(fgets(buf
,STRLEN
,fpin
))
899 if ((temp
=strchr(buf2
,'\n')) != NULL
)
906 if ((temp
=strchr(buf2
,'\n')) != NULL
)
909 if (buf2
[strlen(buf2
)-1]==']')
911 buf2
[strlen(buf2
)-1]='\0';
914 if (gmx_strcasecmp(buf2
,"molecules")==0)
917 fprintf(fpout
,"%s",buf
);
918 } else if (bMolecules
==1)
920 for(i
=0;i
<mtop
->nmolblock
;i
++)
922 nmol
=mtop
->molblock
[i
].nmol
;
923 sprintf(buf
,"%-15s %5d\n",*(mtop
->moltype
[mtop
->molblock
[i
].type
].name
),nmol
);
924 fprintf(fpout
,"%s",buf
);
927 } else if (bMolecules
==2)
932 fprintf(fpout
,"%s",buf
);
936 fprintf(fpout
,"%s",buf
);
941 /* use ffopen to generate backup of topinout */
942 fpout
=ffopen(topfile
,"w");
944 rename(TEMP_FILENM
,topfile
);
948 void md_print_warning(const t_commrec
*cr
,FILE *fplog
,const char *buf
)
952 fprintf(stderr
,"\n%s\n",buf
);
956 fprintf(fplog
,"\n%s\n",buf
);
960 /* simulation conditions to transmit. Keep in mind that they are
961 transmitted to other nodes through an MPI_Reduce after
962 casting them to a real (so the signals can be sent together with other
963 data). This means that the only meaningful values are positive,
965 enum { eglsNABNSB
, eglsCHKPT
, eglsSTOPCOND
, eglsRESETCOUNTERS
, eglsNR
};
966 /* Is the signal in one simulation independent of other simulations? */
967 gmx_bool gs_simlocal
[eglsNR
] = { TRUE
, FALSE
, FALSE
, TRUE
};
970 int nstms
; /* The frequency for intersimulation communication */
971 int sig
[eglsNR
]; /* The signal set by one process in do_md */
972 int set
[eglsNR
]; /* The communicated signal, equal for all processes */
976 static int multisim_min(const gmx_multisim_t
*ms
,int nmin
,int n
)
979 gmx_bool bPos
,bEqual
;
984 gmx_sumi_sim(ms
->nsim
,buf
,ms
);
987 for(s
=0; s
<ms
->nsim
; s
++)
989 bPos
= bPos
&& (buf
[s
] > 0);
990 bEqual
= bEqual
&& (buf
[s
] == buf
[0]);
996 nmin
= min(nmin
,buf
[0]);
1000 /* Find the least common multiple */
1001 for(d
=2; d
<nmin
; d
++)
1004 while (s
< ms
->nsim
&& d
% buf
[s
] == 0)
1010 /* We found the LCM and it is less than nmin */
1022 static int multisim_nstsimsync(const t_commrec
*cr
,
1023 const t_inputrec
*ir
,int repl_ex_nst
)
1030 nmin
= multisim_min(cr
->ms
,nmin
,ir
->nstlist
);
1031 nmin
= multisim_min(cr
->ms
,nmin
,ir
->nstcalcenergy
);
1032 nmin
= multisim_min(cr
->ms
,nmin
,repl_ex_nst
);
1033 if (nmin
== INT_MAX
)
1035 gmx_fatal(FARGS
,"Can not find an appropriate interval for inter-simulation communication, since nstlist, nstcalcenergy and -replex are all <= 0");
1037 /* Avoid inter-simulation communication at every (second) step */
1044 gmx_bcast(sizeof(int),&nmin
,cr
);
1049 static void init_global_signals(globsig_t
*gs
,const t_commrec
*cr
,
1050 const t_inputrec
*ir
,int repl_ex_nst
)
1056 gs
->nstms
= multisim_nstsimsync(cr
,ir
,repl_ex_nst
);
1059 fprintf(debug
,"Syncing simulations for checkpointing and termination every %d steps\n",gs
->nstms
);
1067 for(i
=0; i
<eglsNR
; i
++)
1074 static void copy_coupling_state(t_state
*statea
,t_state
*stateb
,
1075 gmx_ekindata_t
*ekinda
,gmx_ekindata_t
*ekindb
, t_grpopts
* opts
)
1078 /* MRS note -- might be able to get rid of some of the arguments. Look over it when it's all debugged */
1082 /* Make sure we have enough space for x and v */
1083 if (statea
->nalloc
> stateb
->nalloc
)
1085 stateb
->nalloc
= statea
->nalloc
;
1086 srenew(stateb
->x
,stateb
->nalloc
);
1087 srenew(stateb
->v
,stateb
->nalloc
);
1090 stateb
->natoms
= statea
->natoms
;
1091 stateb
->ngtc
= statea
->ngtc
;
1092 stateb
->nnhpres
= statea
->nnhpres
;
1093 stateb
->veta
= statea
->veta
;
1096 copy_mat(ekinda
->ekin
,ekindb
->ekin
);
1097 for (i
=0; i
<stateb
->ngtc
; i
++)
1099 ekindb
->tcstat
[i
].T
= ekinda
->tcstat
[i
].T
;
1100 ekindb
->tcstat
[i
].Th
= ekinda
->tcstat
[i
].Th
;
1101 copy_mat(ekinda
->tcstat
[i
].ekinh
,ekindb
->tcstat
[i
].ekinh
);
1102 copy_mat(ekinda
->tcstat
[i
].ekinf
,ekindb
->tcstat
[i
].ekinf
);
1103 ekindb
->tcstat
[i
].ekinscalef_nhc
= ekinda
->tcstat
[i
].ekinscalef_nhc
;
1104 ekindb
->tcstat
[i
].ekinscaleh_nhc
= ekinda
->tcstat
[i
].ekinscaleh_nhc
;
1105 ekindb
->tcstat
[i
].vscale_nhc
= ekinda
->tcstat
[i
].vscale_nhc
;
1108 copy_rvecn(statea
->x
,stateb
->x
,0,stateb
->natoms
);
1109 copy_rvecn(statea
->v
,stateb
->v
,0,stateb
->natoms
);
1110 copy_mat(statea
->box
,stateb
->box
);
1111 copy_mat(statea
->box_rel
,stateb
->box_rel
);
1112 copy_mat(statea
->boxv
,stateb
->boxv
);
1114 for (i
= 0; i
<stateb
->ngtc
; i
++)
1116 nc
= i
*opts
->nhchainlength
;
1117 for (j
=0; j
<opts
->nhchainlength
; j
++)
1119 stateb
->nosehoover_xi
[nc
+j
] = statea
->nosehoover_xi
[nc
+j
];
1120 stateb
->nosehoover_vxi
[nc
+j
] = statea
->nosehoover_vxi
[nc
+j
];
1123 if (stateb
->nhpres_xi
!= NULL
)
1125 for (i
= 0; i
<stateb
->nnhpres
; i
++)
1127 nc
= i
*opts
->nhchainlength
;
1128 for (j
=0; j
<opts
->nhchainlength
; j
++)
1130 stateb
->nhpres_xi
[nc
+j
] = statea
->nhpres_xi
[nc
+j
];
1131 stateb
->nhpres_vxi
[nc
+j
] = statea
->nhpres_vxi
[nc
+j
];
1137 static void compute_globals(FILE *fplog
, gmx_global_stat_t gstat
, t_commrec
*cr
, t_inputrec
*ir
,
1138 t_forcerec
*fr
, gmx_ekindata_t
*ekind
,
1139 t_state
*state
, t_state
*state_global
, t_mdatoms
*mdatoms
,
1140 t_nrnb
*nrnb
, t_vcm
*vcm
, gmx_wallcycle_t wcycle
,
1141 gmx_enerdata_t
*enerd
,tensor force_vir
, tensor shake_vir
, tensor total_vir
,
1142 tensor pres
, rvec mu_tot
, gmx_constr_t constr
,
1143 globsig_t
*gs
,gmx_bool bInterSimGS
,
1144 matrix box
, gmx_mtop_t
*top_global
, real
*pcurr
,
1145 int natoms
, gmx_bool
*bSumEkinhOld
, int flags
)
1148 real gs_buf
[eglsNR
];
1149 tensor corr_vir
,corr_pres
;
1150 gmx_bool bEner
,bPres
,bTemp
;
1151 gmx_bool bRerunMD
, bStopCM
, bGStat
, bIterate
,
1152 bFirstIterate
,bReadEkin
,bEkinAveVel
,bScaleEkin
, bConstrain
;
1153 real prescorr
,enercorr
,dvdlcorr
;
1155 /* translate CGLO flags to gmx_booleans */
1156 bRerunMD
= flags
& CGLO_RERUNMD
;
1157 bStopCM
= flags
& CGLO_STOPCM
;
1158 bGStat
= flags
& CGLO_GSTAT
;
1159 /* FIX ME after 4.5 */
1160 /* temporary hack because we are using gmx_bool (unsigned char) */
1162 bReadEkin
= (flags
& CGLO_READEKIN
) != 0;
1163 bScaleEkin
= (flags
& CGLO_SCALEEKIN
) != 0;
1164 bEner
= flags
& CGLO_ENERGY
;
1165 bTemp
= flags
& CGLO_TEMPERATURE
;
1166 bPres
= (flags
& CGLO_PRESSURE
) != 0;
1167 bConstrain
= (flags
& CGLO_CONSTRAINT
) != 0;
1168 bIterate
= (flags
& CGLO_ITERATE
) != 0;
1169 bFirstIterate
= (flags
& CGLO_FIRSTITERATE
) != 0;
1171 /* we calculate a full state kinetic energy either with full-step velocity verlet
1172 or half step where we need the pressure */
1173 bEkinAveVel
= (ir
->eI
==eiVV
|| (ir
->eI
==eiVVAK
&& IR_NPT_TROTTER(ir
) && bPres
) || bReadEkin
);
1175 /* in initalization, it sums the shake virial in vv, and to
1176 sums ekinh_old in leapfrog (or if we are calculating ekinh_old for other reasons */
1178 /* ########## Kinetic energy ############## */
1182 /* Non-equilibrium MD: this is parallellized, but only does communication
1183 * when there really is NEMD.
1186 if (PAR(cr
) && (ekind
->bNEMD
))
1188 accumulate_u(cr
,&(ir
->opts
),ekind
);
1193 restore_ekinstate_from_state(cr
,ekind
,&state_global
->ekinstate
);
1198 calc_ke_part(state
,&(ir
->opts
),mdatoms
,ekind
,nrnb
,bEkinAveVel
,bIterate
);
1203 /* Calculate center of mass velocity if necessary, also parallellized */
1204 if (bStopCM
&& !bRerunMD
&& bEner
)
1206 calc_vcm_grp(fplog
,mdatoms
->start
,mdatoms
->homenr
,mdatoms
,
1207 state
->x
,state
->v
,vcm
);
1211 if (bTemp
|| bPres
|| bEner
|| bConstrain
)
1215 /* We will not sum ekinh_old,
1216 * so signal that we still have to do it.
1218 *bSumEkinhOld
= TRUE
;
1225 for(i
=0; i
<eglsNR
; i
++)
1227 gs_buf
[i
] = gs
->sig
[i
];
1232 wallcycle_start(wcycle
,ewcMoveE
);
1233 GMX_MPE_LOG(ev_global_stat_start
);
1234 global_stat(fplog
,gstat
,cr
,enerd
,force_vir
,shake_vir
,mu_tot
,
1235 ir
,ekind
,constr
,vcm
,
1236 gs
!= NULL
? eglsNR
: 0,gs_buf
,
1238 *bSumEkinhOld
,flags
);
1239 GMX_MPE_LOG(ev_global_stat_finish
);
1240 wallcycle_stop(wcycle
,ewcMoveE
);
1244 if (MULTISIM(cr
) && bInterSimGS
)
1248 /* Communicate the signals between the simulations */
1249 gmx_sum_sim(eglsNR
,gs_buf
,cr
->ms
);
1251 /* Communicate the signals form the master to the others */
1252 gmx_bcast(eglsNR
*sizeof(gs_buf
[0]),gs_buf
,cr
);
1254 for(i
=0; i
<eglsNR
; i
++)
1256 if (bInterSimGS
|| gs_simlocal
[i
])
1258 /* Set the communicated signal only when it is non-zero,
1259 * since signals might not be processed at each MD step.
1261 gsi
= (gs_buf
[i
] >= 0 ?
1262 (int)(gs_buf
[i
] + 0.5) :
1263 (int)(gs_buf
[i
] - 0.5));
1268 /* Turn off the local signal */
1273 *bSumEkinhOld
= FALSE
;
1277 if (!ekind
->bNEMD
&& debug
&& bTemp
&& (vcm
->nr
> 0))
1280 mdatoms
->start
,mdatoms
->start
+mdatoms
->homenr
,
1281 state
->v
,vcm
->group_p
[0],
1282 mdatoms
->massT
,mdatoms
->tmass
,ekind
->ekin
);
1286 /* Do center of mass motion removal */
1287 if (bStopCM
&& !bRerunMD
) /* is this correct? Does it get called too often with this logic? */
1289 check_cm_grp(fplog
,vcm
,ir
,1);
1290 do_stopcm_grp(fplog
,mdatoms
->start
,mdatoms
->homenr
,mdatoms
->cVCM
,
1291 state
->x
,state
->v
,vcm
);
1292 inc_nrnb(nrnb
,eNR_STOPCM
,mdatoms
->homenr
);
1298 /* Sum the kinetic energies of the groups & calc temp */
1299 /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with IR_NPT_TROTTER */
1300 /* three maincase: VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
1301 Leap with AveVel is also an option for the future but not supported now.
1302 bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
1303 If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
1304 bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
1305 If FALSE, we go ahead and erase over it.
1307 enerd
->term
[F_TEMP
] = sum_ekin(&(ir
->opts
),ekind
,&(enerd
->term
[F_DKDL
]),
1308 bEkinAveVel
,bIterate
,bScaleEkin
);
1310 enerd
->term
[F_EKIN
] = trace(ekind
->ekin
);
1313 /* ########## Long range energy information ###### */
1315 if (bEner
|| bPres
|| bConstrain
)
1317 calc_dispcorr(fplog
,ir
,fr
,0,top_global
->natoms
,box
,state
->lambda
,
1318 corr_pres
,corr_vir
,&prescorr
,&enercorr
,&dvdlcorr
);
1321 if (bEner
&& bFirstIterate
)
1323 enerd
->term
[F_DISPCORR
] = enercorr
;
1324 enerd
->term
[F_EPOT
] += enercorr
;
1325 enerd
->term
[F_DVDL
] += dvdlcorr
;
1326 if (fr
->efep
!= efepNO
) {
1327 enerd
->dvdl_lin
+= dvdlcorr
;
1331 /* ########## Now pressure ############## */
1332 if (bPres
|| bConstrain
)
1335 m_add(force_vir
,shake_vir
,total_vir
);
1337 /* Calculate pressure and apply LR correction if PPPM is used.
1338 * Use the box from last timestep since we already called update().
1341 enerd
->term
[F_PRES
] = calc_pres(fr
->ePBC
,ir
->nwall
,box
,ekind
->ekin
,total_vir
,pres
,
1342 (fr
->eeltype
==eelPPPM
)?enerd
->term
[F_COUL_RECIP
]:0.0);
1344 /* Calculate long range corrections to pressure and energy */
1345 /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
1346 and computes enerd->term[F_DISPCORR]. Also modifies the
1347 total_vir and pres tesors */
1349 m_add(total_vir
,corr_vir
,total_vir
);
1350 m_add(pres
,corr_pres
,pres
);
1351 enerd
->term
[F_PDISPCORR
] = prescorr
;
1352 enerd
->term
[F_PRES
] += prescorr
;
1353 *pcurr
= enerd
->term
[F_PRES
];
1354 /* calculate temperature using virial */
1355 enerd
->term
[F_VTEMP
] = calc_temp(trace(total_vir
),ir
->opts
.nrdf
[0]);
1361 /* Definitions for convergence of iterated constraints */
1363 /* iterate constraints up to 50 times */
1364 #define MAXITERCONST 50
1369 real f
,fprev
,x
,xprev
;
1372 real allrelerr
[MAXITERCONST
+2];
1373 int num_close
; /* number of "close" violations, caused by limited precision. */
1377 #define CONVERGEITER 0.000000001
1378 #define CLOSE_ENOUGH 0.000001000
1380 #define CONVERGEITER 0.0001
1381 #define CLOSE_ENOUGH 0.0050
1384 /* we want to keep track of the close calls. If there are too many, there might be some other issues.
1385 so we make sure that it's either less than some predetermined number, or if more than that number,
1386 only some small fraction of the total. */
1387 #define MAX_NUMBER_CLOSE 50
1388 #define FRACTION_CLOSE 0.001
1390 /* maximum length of cyclic traps to check, emerging from limited numerical precision */
1393 static void gmx_iterate_init(gmx_iterate_t
*iterate
,gmx_bool bIterate
)
1397 iterate
->iter_i
= 0;
1398 iterate
->bIterate
= bIterate
;
1399 iterate
->num_close
= 0;
1400 for (i
=0;i
<MAXITERCONST
+2;i
++)
1402 iterate
->allrelerr
[i
] = 0;
1406 static gmx_bool
done_iterating(const t_commrec
*cr
,FILE *fplog
, int nsteps
, gmx_iterate_t
*iterate
, gmx_bool bFirstIterate
, real fom
, real
*newf
)
1408 /* monitor convergence, and use a secant search to propose new
1411 The secant method computes x_{i+1} = x_{i} - f(x_{i}) * ---------------------
1412 f(x_{i}) - f(x_{i-1})
1414 The function we are trying to zero is fom-x, where fom is the
1415 "figure of merit" which is the pressure (or the veta value) we
1416 would get by putting in an old value of the pressure or veta into
1417 the incrementor function for the step or half step. I have
1418 verified that this gives the same answer as self consistent
1419 iteration, usually in many fewer steps, especially for small tau_p.
1421 We could possibly eliminate an iteration with proper use
1422 of the value from the previous step, but that would take a bit
1423 more bookkeeping, especially for veta, since tests indicate the
1424 function of veta on the last step is not sufficiently close to
1425 guarantee convergence this step. This is
1426 good enough for now. On my tests, I could use tau_p down to
1427 0.02, which is smaller that would ever be necessary in
1428 practice. Generally, 3-5 iterations will be sufficient */
1438 iterate
->f
= fom
-iterate
->x
;
1445 iterate
->f
= fom
-iterate
->x
; /* we want to zero this difference */
1446 if ((iterate
->iter_i
> 1) && (iterate
->iter_i
< MAXITERCONST
))
1448 if (iterate
->f
==iterate
->fprev
)
1454 *newf
= iterate
->x
- (iterate
->x
-iterate
->xprev
)*(iterate
->f
)/(iterate
->f
-iterate
->fprev
);
1459 /* just use self-consistent iteration the first step to initialize, or
1460 if it's not converging (which happens occasionally -- need to investigate why) */
1464 /* Consider a slight shortcut allowing us to exit one sooner -- we check the
1465 difference between the closest of x and xprev to the new
1466 value. To be 100% certain, we should check the difference between
1467 the last result, and the previous result, or
1469 relerr = (fabs((x-xprev)/fom));
1471 but this is pretty much never necessary under typical conditions.
1472 Checking numerically, it seems to lead to almost exactly the same
1473 trajectories, but there are small differences out a few decimal
1474 places in the pressure, and eventually in the v_eta, but it could
1477 if (fabs(*newf-x) < fabs(*newf - xprev)) { xmin = x;} else { xmin = xprev;}
1478 relerr = (fabs((*newf-xmin) / *newf));
1481 err
= fabs((iterate
->f
-iterate
->fprev
));
1482 relerr
= fabs(err
/fom
);
1484 iterate
->allrelerr
[iterate
->iter_i
] = relerr
;
1486 if (iterate
->iter_i
> 0)
1490 fprintf(debug
,"Iterating NPT constraints: %6i %20.12f%14.6g%20.12f\n",
1491 iterate
->iter_i
,fom
,relerr
,*newf
);
1494 if ((relerr
< CONVERGEITER
) || (err
< CONVERGEITER
) || (fom
==0) || ((iterate
->x
== iterate
->xprev
) && iterate
->iter_i
> 1))
1496 iterate
->bIterate
= FALSE
;
1499 fprintf(debug
,"Iterating NPT constraints: CONVERGED\n");
1503 if (iterate
->iter_i
> MAXITERCONST
)
1505 if (relerr
< CLOSE_ENOUGH
)
1508 for (i
=1;i
<CYCLEMAX
;i
++) {
1509 if ((iterate
->allrelerr
[iterate
->iter_i
-(1+i
)] == iterate
->allrelerr
[iterate
->iter_i
-1]) &&
1510 (iterate
->allrelerr
[iterate
->iter_i
-(1+i
)] == iterate
->allrelerr
[iterate
->iter_i
-(1+2*i
)])) {
1514 fprintf(debug
,"Exiting from an NPT iterating cycle of length %d\n",i
);
1521 /* step 1: trapped in a numerical attractor */
1522 /* we are trapped in a numerical attractor, and can't converge any more, and are close to the final result.
1523 Better to give up convergence here than have the simulation die.
1525 iterate
->num_close
++;
1530 /* Step #2: test if we are reasonably close for other reasons, then monitor the number. If not, die */
1532 /* how many close calls have we had? If less than a few, we're OK */
1533 if (iterate
->num_close
< MAX_NUMBER_CLOSE
)
1535 sprintf(buf
,"Slight numerical convergence deviation with NPT at step %d, relative error only %10.5g, likely not a problem, continuing\n",nsteps
,relerr
);
1536 md_print_warning(cr
,fplog
,buf
);
1537 iterate
->num_close
++;
1539 /* if more than a few, check the total fraction. If too high, die. */
1540 } else if (iterate
->num_close
/(double)nsteps
> FRACTION_CLOSE
) {
1541 gmx_fatal(FARGS
,"Could not converge NPT constraints, too many exceptions (%d%%\n",iterate
->num_close
/(double)nsteps
);
1547 gmx_fatal(FARGS
,"Could not converge NPT constraints\n");
1552 iterate
->xprev
= iterate
->x
;
1554 iterate
->fprev
= iterate
->f
;
1560 static void check_nst_param(FILE *fplog
,t_commrec
*cr
,
1561 const char *desc_nst
,int nst
,
1562 const char *desc_p
,int *p
)
1566 if (*p
> 0 && *p
% nst
!= 0)
1568 /* Round up to the next multiple of nst */
1569 *p
= ((*p
)/nst
+ 1)*nst
;
1570 sprintf(buf
,"NOTE: %s changes %s to %d\n",desc_nst
,desc_p
,*p
);
1571 md_print_warning(cr
,fplog
,buf
);
1575 static void reset_all_counters(FILE *fplog
,t_commrec
*cr
,
1576 gmx_large_int_t step
,
1577 gmx_large_int_t
*step_rel
,t_inputrec
*ir
,
1578 gmx_wallcycle_t wcycle
,t_nrnb
*nrnb
,
1579 gmx_runtime_t
*runtime
)
1581 char buf
[STRLEN
],sbuf
[STEPSTRSIZE
];
1583 /* Reset all the counters related to performance over the run */
1584 sprintf(buf
,"Step %s: resetting all time and cycle counters\n",
1585 gmx_step_str(step
,sbuf
));
1586 md_print_warning(cr
,fplog
,buf
);
1588 wallcycle_stop(wcycle
,ewcRUN
);
1589 wallcycle_reset_all(wcycle
);
1590 if (DOMAINDECOMP(cr
))
1592 reset_dd_statistics_counters(cr
->dd
);
1595 ir
->init_step
+= *step_rel
;
1596 ir
->nsteps
-= *step_rel
;
1598 wallcycle_start(wcycle
,ewcRUN
);
1599 runtime_start(runtime
);
1600 print_date_and_time(fplog
,cr
->nodeid
,"Restarted time",runtime
);
1603 static int check_nstglobalcomm(FILE *fplog
,t_commrec
*cr
,
1604 int nstglobalcomm
,t_inputrec
*ir
)
1608 if (!EI_DYNAMICS(ir
->eI
))
1613 if (nstglobalcomm
== -1)
1615 if (ir
->nstcalcenergy
== 0 && ir
->nstlist
== 0)
1618 if (ir
->nstenergy
> 0 && ir
->nstenergy
< nstglobalcomm
)
1620 nstglobalcomm
= ir
->nstenergy
;
1625 /* We assume that if nstcalcenergy > nstlist,
1626 * nstcalcenergy is a multiple of nstlist.
1628 if (ir
->nstcalcenergy
== 0 ||
1629 (ir
->nstlist
> 0 && ir
->nstlist
< ir
->nstcalcenergy
))
1631 nstglobalcomm
= ir
->nstlist
;
1635 nstglobalcomm
= ir
->nstcalcenergy
;
1641 if (ir
->nstlist
> 0 &&
1642 nstglobalcomm
> ir
->nstlist
&& nstglobalcomm
% ir
->nstlist
!= 0)
1644 nstglobalcomm
= (nstglobalcomm
/ ir
->nstlist
)*ir
->nstlist
;
1645 sprintf(buf
,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm
);
1646 md_print_warning(cr
,fplog
,buf
);
1648 if (nstglobalcomm
> ir
->nstcalcenergy
)
1650 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
1651 "nstcalcenergy",&ir
->nstcalcenergy
);
1654 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
1655 "nstenergy",&ir
->nstenergy
);
1657 check_nst_param(fplog
,cr
,"-gcom",nstglobalcomm
,
1658 "nstlog",&ir
->nstlog
);
1661 if (ir
->comm_mode
!= ecmNO
&& ir
->nstcomm
< nstglobalcomm
)
1663 sprintf(buf
,"WARNING: Changing nstcomm from %d to %d\n",
1664 ir
->nstcomm
,nstglobalcomm
);
1665 md_print_warning(cr
,fplog
,buf
);
1666 ir
->nstcomm
= nstglobalcomm
;
1669 return nstglobalcomm
;
1672 void check_ir_old_tpx_versions(t_commrec
*cr
,FILE *fplog
,
1673 t_inputrec
*ir
,gmx_mtop_t
*mtop
)
1675 /* Check required for old tpx files */
1676 if (IR_TWINRANGE(*ir
) && ir
->nstlist
> 1 &&
1677 ir
->nstcalcenergy
% ir
->nstlist
!= 0)
1679 md_print_warning(cr
,fplog
,"Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies");
1681 if (gmx_mtop_ftype_count(mtop
,F_CONSTR
) +
1682 gmx_mtop_ftype_count(mtop
,F_CONSTRNC
) > 0 &&
1683 ir
->eConstrAlg
== econtSHAKE
)
1685 md_print_warning(cr
,fplog
,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect");
1686 if (ir
->epc
!= epcNO
)
1688 gmx_fatal(FARGS
,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
1691 check_nst_param(fplog
,cr
,"nstlist",ir
->nstlist
,
1692 "nstcalcenergy",&ir
->nstcalcenergy
);
1693 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
1694 "nstenergy",&ir
->nstenergy
);
1695 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
1696 "nstlog",&ir
->nstlog
);
1697 if (ir
->efep
!= efepNO
)
1699 check_nst_param(fplog
,cr
,"nstcalcenergy",ir
->nstcalcenergy
,
1700 "nstdhdl",&ir
->nstdhdl
);
1706 gmx_bool bGStatEveryStep
;
1707 gmx_large_int_t step_ns
;
1708 gmx_large_int_t step_nscheck
;
1709 gmx_large_int_t nns
;
1719 static void reset_nlistheuristics(gmx_nlheur_t
*nlh
,gmx_large_int_t step
)
1723 nlh
->step_nscheck
= step
;
1726 static void init_nlistheuristics(gmx_nlheur_t
*nlh
,
1727 gmx_bool bGStatEveryStep
,gmx_large_int_t step
)
1729 nlh
->bGStatEveryStep
= bGStatEveryStep
;
1736 reset_nlistheuristics(nlh
,step
);
1739 static void update_nliststatistics(gmx_nlheur_t
*nlh
,gmx_large_int_t step
)
1741 gmx_large_int_t nl_lt
;
1742 char sbuf
[STEPSTRSIZE
],sbuf2
[STEPSTRSIZE
];
1744 /* Determine the neighbor list life time */
1745 nl_lt
= step
- nlh
->step_ns
;
1748 fprintf(debug
,"%d atoms beyond ns buffer, updating neighbor list after %s steps\n",nlh
->nabnsb
,gmx_step_str(nl_lt
,sbuf
));
1752 nlh
->s2
+= nl_lt
*nl_lt
;
1753 nlh
->ab
+= nlh
->nabnsb
;
1754 if (nlh
->lt_runav
== 0)
1756 nlh
->lt_runav
= nl_lt
;
1757 /* Initialize the fluctuation average
1758 * such that at startup we check after 0 steps.
1760 nlh
->lt_runav2
= sqr(nl_lt
/2.0);
1762 /* Running average with 0.9 gives an exp. history of 9.5 */
1763 nlh
->lt_runav2
= 0.9*nlh
->lt_runav2
+ 0.1*sqr(nlh
->lt_runav
- nl_lt
);
1764 nlh
->lt_runav
= 0.9*nlh
->lt_runav
+ 0.1*nl_lt
;
1765 if (nlh
->bGStatEveryStep
)
1767 /* Always check the nlist validity */
1768 nlh
->step_nscheck
= step
;
1772 /* We check after: <life time> - 2*sigma
1773 * The factor 2 is quite conservative,
1774 * but we assume that with nstlist=-1 the user
1775 * prefers exact integration over performance.
1777 nlh
->step_nscheck
= step
1778 + (int)(nlh
->lt_runav
- 2.0*sqrt(nlh
->lt_runav2
)) - 1;
1782 fprintf(debug
,"nlist life time %s run av. %4.1f sig %3.1f check %s check with -gcom %d\n",
1783 gmx_step_str(nl_lt
,sbuf
),nlh
->lt_runav
,sqrt(nlh
->lt_runav2
),
1784 gmx_step_str(nlh
->step_nscheck
-step
+1,sbuf2
),
1785 (int)(nlh
->lt_runav
- 2.0*sqrt(nlh
->lt_runav2
)));
1789 static void set_nlistheuristics(gmx_nlheur_t
*nlh
,gmx_bool bReset
,gmx_large_int_t step
)
1795 reset_nlistheuristics(nlh
,step
);
1799 update_nliststatistics(nlh
,step
);
1802 nlh
->step_ns
= step
;
1803 /* Initialize the cumulative coordinate scaling matrix */
1804 clear_mat(nlh
->scale_tot
);
1805 for(d
=0; d
<DIM
; d
++)
1807 nlh
->scale_tot
[d
][d
] = 1.0;
1811 double do_md_membed(FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
1812 const output_env_t oenv
, gmx_bool bVerbose
,gmx_bool bCompact
,
1814 gmx_vsite_t
*vsite
,gmx_constr_t constr
,
1815 int stepout
,t_inputrec
*ir
,
1816 gmx_mtop_t
*top_global
,
1818 t_state
*state_global
,
1820 t_nrnb
*nrnb
,gmx_wallcycle_t wcycle
,
1821 gmx_edsam_t ed
,t_forcerec
*fr
,
1822 int repl_ex_nst
,int repl_ex_seed
,
1823 real cpt_period
,real max_hours
,
1824 const char *deviceOptions
,
1825 unsigned long Flags
,
1826 gmx_runtime_t
*runtime
,
1827 rvec fac
, rvec
*r_ins
, pos_ins_t
*pos_ins
, t_block
*ins_at
,
1828 real xy_step
, real z_step
, int it_xy
, int it_z
)
1831 gmx_large_int_t step
,step_rel
;
1834 gmx_bool bGStatEveryStep
,bGStat
,bNstEner
,bCalcEnerPres
;
1835 gmx_bool bNS
,bNStList
,bSimAnn
,bStopCM
,bRerunMD
,bNotLastFrame
=FALSE
,
1836 bFirstStep
,bStateFromTPX
,bInitStep
,bLastStep
,
1837 bBornRadii
,bStartingFromCpt
;
1838 gmx_bool bDoDHDL
=FALSE
;
1839 gmx_bool do_ene
,do_log
,do_verbose
,bRerunWarnNoV
=TRUE
,
1840 bForceUpdate
=FALSE
,bCPT
;
1842 gmx_bool bMasterState
;
1843 int force_flags
,cglo_flags
;
1844 tensor force_vir
,shake_vir
,total_vir
,tmp_vir
,pres
;
1846 t_trxstatus
*status
;
1849 t_state
*bufstate
=NULL
;
1850 matrix
*scale_tot
,pcoupl_mu
,M
,ebox
;
1852 t_trxframe rerun_fr
;
1853 /* gmx_repl_ex_t repl_ex=NULL;*/
1856 gmx_localtop_t
*top
;
1857 t_mdebin
*mdebin
=NULL
;
1858 t_state
*state
=NULL
;
1859 rvec
*f_global
=NULL
;
1862 gmx_enerdata_t
*enerd
;
1864 gmx_global_stat_t gstat
;
1865 gmx_update_t upd
=NULL
;
1866 t_graph
*graph
=NULL
;
1870 gmx_groups_t
*groups
;
1871 gmx_ekindata_t
*ekind
, *ekind_save
;
1872 gmx_shellfc_t shellfc
;
1873 int count
,nconverged
=0;
1876 gmx_bool bIonize
=FALSE
;
1877 gmx_bool bTCR
=FALSE
,bConverged
=TRUE
,bOK
,bSumEkinhOld
,bExchanged
;
1879 gmx_bool bResetCountersHalfMaxH
=FALSE
;
1880 gmx_bool bVV
,bIterations
,bIterate
,bFirstIterate
,bTemp
,bPres
,bTrotter
;
1883 rvec
*xcopy
=NULL
,*vcopy
=NULL
,*cbuf
=NULL
;
1884 matrix boxcopy
={{0}},lastbox
;
1885 real veta_save
,pcurr
,scalevir
,tracevir
;
1888 real last_conserved
= 0;
1892 char sbuf
[STEPSTRSIZE
],sbuf2
[STEPSTRSIZE
];
1893 int handled_stop_condition
=gmx_stop_cond_none
; /* compare to get_stop_condition*/
1894 gmx_iterate_t iterate
;
1896 /* Temporary addition for FAHCORE checkpointing */
1900 /* Check for special mdrun options */
1901 bRerunMD
= (Flags
& MD_RERUN
);
1902 bIonize
= (Flags
& MD_IONIZE
);
1903 bFFscan
= (Flags
& MD_FFSCAN
);
1904 bAppend
= (Flags
& MD_APPENDFILES
);
1905 bGStatEveryStep
= FALSE
;
1906 if (Flags
& MD_RESETCOUNTERSHALFWAY
)
1910 /* Signal to reset the counters half the simulation steps. */
1911 wcycle_set_reset_counters(wcycle
,ir
->nsteps
/2);
1913 /* Signal to reset the counters halfway the simulation time. */
1914 bResetCountersHalfMaxH
= (max_hours
> 0);
1917 /* md-vv uses averaged full step velocities for T-control
1918 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
1919 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
1920 bVV
= EI_VV(ir
->eI
);
1921 if (bVV
) /* to store the initial velocities while computing virial */
1923 snew(cbuf
,top_global
->natoms
);
1925 /* all the iteratative cases - only if there are constraints */
1926 bIterations
= ((IR_NPT_TROTTER(ir
)) && (constr
) && (!bRerunMD
));
1927 bTrotter
= (bVV
&& (IR_NPT_TROTTER(ir
) || (IR_NVT_TROTTER(ir
))));
1931 /* Since we don't know if the frames read are related in any way,
1932 * rebuild the neighborlist at every step.
1935 ir
->nstcalcenergy
= 1;
1939 check_ir_old_tpx_versions(cr
,fplog
,ir
,top_global
);
1941 nstglobalcomm
= check_nstglobalcomm(fplog
,cr
,nstglobalcomm
,ir
);
1942 /*bGStatEveryStep = (nstglobalcomm == 1);*/
1943 bGStatEveryStep
= FALSE
;
1945 if (!bGStatEveryStep
&& ir
->nstlist
== -1 && fplog
!= NULL
)
1948 "To reduce the energy communication with nstlist = -1\n"
1949 "the neighbor list validity should not be checked at every step,\n"
1950 "this means that exact integration is not guaranteed.\n"
1951 "The neighbor list validity is checked after:\n"
1952 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
1953 "In most cases this will result in exact integration.\n"
1954 "This reduces the energy communication by a factor of 2 to 3.\n"
1955 "If you want less energy communication, set nstlist > 3.\n\n");
1958 if (bRerunMD
|| bFFscan
)
1962 groups
= &top_global
->groups
;
1964 /* Initial values */
1965 init_md(fplog
,cr
,ir
,oenv
,&t
,&t0
,&state_global
->lambda
,&lam0
,
1966 nrnb
,top_global
,&upd
,
1967 nfile
,fnm
,&outf
,&mdebin
,
1968 force_vir
,shake_vir
,mu_tot
,&bSimAnn
,&vcm
,state_global
,Flags
);
1970 clear_mat(total_vir
);
1972 /* Energy terms and groups */
1974 init_enerdata(top_global
->groups
.grps
[egcENER
].nr
,ir
->n_flambda
,enerd
);
1975 if (DOMAINDECOMP(cr
))
1981 snew(f
,top_global
->natoms
);
1984 /* Kinetic energy data */
1986 init_ekindata(fplog
,top_global
,&(ir
->opts
),ekind
);
1987 /* needed for iteration of constraints */
1989 init_ekindata(fplog
,top_global
,&(ir
->opts
),ekind_save
);
1990 /* Copy the cos acceleration to the groups struct */
1991 ekind
->cosacc
.cos_accel
= ir
->cos_accel
;
1993 gstat
= global_stat_init(ir
);
1996 /* Check for polarizable models and flexible constraints */
1997 shellfc
= init_shell_flexcon(fplog
,
1998 top_global
,n_flexible_constraints(constr
),
1999 (ir
->bContinuation
||
2000 (DOMAINDECOMP(cr
) && !MASTER(cr
))) ?
2001 NULL
: state_global
->x
);
2006 tMPI_Thread_mutex_lock(&deform_init_box_mutex);
2008 set_deform_reference_box(upd,
2009 deform_init_init_step_tpx,
2010 deform_init_box_tpx);
2012 tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
2017 double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
2018 if ((io > 2000) && MASTER(cr))
2020 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
2024 if (DOMAINDECOMP(cr
)) {
2025 top
= dd_init_local_top(top_global
);
2028 dd_init_local_state(cr
->dd
,state_global
,state
);
2030 if (DDMASTER(cr
->dd
) && ir
->nstfout
) {
2031 snew(f_global
,state_global
->natoms
);
2035 /* Initialize the particle decomposition and split the topology */
2036 top
= split_system(fplog
,top_global
,ir
,cr
);
2038 pd_cg_range(cr
,&fr
->cg0
,&fr
->hcg
);
2039 pd_at_range(cr
,&a0
,&a1
);
2041 top
= gmx_mtop_generate_local_top(top_global
,ir
);
2044 a1
= top_global
->natoms
;
2047 state
= partdec_init_local_state(cr
,state_global
);
2050 atoms2md(top_global
,ir
,0,NULL
,a0
,a1
-a0
,mdatoms
);
2053 set_vsite_top(vsite
,top
,mdatoms
,cr
);
2056 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
) {
2057 graph
= mk_graph(fplog
,&(top
->idef
),0,top_global
->natoms
,FALSE
,FALSE
);
2061 make_local_shells(cr
,mdatoms
,shellfc
);
2064 if (ir
->pull
&& PAR(cr
)) {
2065 dd_make_local_pull_groups(NULL
,ir
->pull
,mdatoms
);
2069 if (DOMAINDECOMP(cr
))
2071 /* Distribute the charge groups over the nodes from the master node */
2072 dd_partition_system(fplog
,ir
->init_step
,cr
,TRUE
,1,
2073 state_global
,top_global
,ir
,
2074 state
,&f
,mdatoms
,top
,fr
,
2075 vsite
,shellfc
,constr
,
2079 update_mdatoms(mdatoms
,state
->lambda
);
2083 /* Update mdebin with energy history if appending to output files */
2084 if ( Flags
& MD_APPENDFILES
)
2086 restore_energyhistory_from_state(mdebin
,&state_global
->enerhist
);
2088 /* Set the initial energy history in state to zero by updating once */
2089 update_energyhistory(&state_global
->enerhist
,mdebin
);
2092 if ((state
->flags
& (1<<estLD_RNG
)) && (Flags
& MD_READ_RNG
)) {
2093 /* Set the random state if we read a checkpoint file */
2094 set_stochd_state(upd
,state
);
2097 /* Initialize constraints */
2099 if (!DOMAINDECOMP(cr
))
2100 set_constraints(constr
,top
,ir
,mdatoms
,cr
);
2103 /* Check whether we have to GCT stuff */
2104 /* bTCR = ftp2bSet(efGCT,nfile,fnm);
2107 fprintf(stderr,"Will do General Coupling Theory!\n");
2109 gnx = top_global->mols.nr;
2111 for(i=0; (i<gnx); i++) {
2116 /* if (repl_ex_nst > 0 && MASTER(cr))
2117 repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
2118 repl_ex_nst,repl_ex_seed);*/
2120 if (!ir
->bContinuation
&& !bRerunMD
)
2122 if (mdatoms
->cFREEZE
&& (state
->flags
& (1<<estV
)))
2124 /* Set the velocities of frozen particles to zero */
2125 for(i
=mdatoms
->start
; i
<mdatoms
->start
+mdatoms
->homenr
; i
++)
2127 for(m
=0; m
<DIM
; m
++)
2129 if (ir
->opts
.nFreeze
[mdatoms
->cFREEZE
[i
]][m
])
2139 /* Constrain the initial coordinates and velocities */
2140 do_constrain_first(fplog
,constr
,ir
,mdatoms
,state
,f
,
2141 graph
,cr
,nrnb
,fr
,top
,shake_vir
);
2145 /* Construct the virtual sites for the initial configuration */
2146 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,NULL
,
2147 top
->idef
.iparams
,top
->idef
.il
,
2148 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
2154 /* I'm assuming we need global communication the first time! MRS */
2155 cglo_flags
= (CGLO_TEMPERATURE
| CGLO_GSTAT
2156 | (bVV
? CGLO_PRESSURE
:0)
2157 | (bVV
? CGLO_CONSTRAINT
:0)
2158 | (bRerunMD
? CGLO_RERUNMD
:0)
2159 | ((Flags
& MD_READ_EKIN
) ? CGLO_READEKIN
:0));
2161 bSumEkinhOld
= FALSE
;
2162 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2163 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
2164 constr
,NULL
,FALSE
,state
->box
,
2165 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,cglo_flags
);
2166 if (ir
->eI
== eiVVAK
) {
2167 /* a second call to get the half step temperature initialized as well */
2168 /* we do the same call as above, but turn the pressure off -- internally, this
2169 is recognized as a velocity verlet half-step kinetic energy calculation.
2170 This minimized excess variables, but perhaps loses some logic?*/
2172 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2173 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
2174 constr
,NULL
,FALSE
,state
->box
,
2175 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2176 cglo_flags
&~ CGLO_PRESSURE
);
2179 /* Calculate the initial half step temperature, and save the ekinh_old */
2180 if (!(Flags
& MD_STARTFROMCPT
))
2182 for(i
=0; (i
<ir
->opts
.ngtc
); i
++)
2184 copy_mat(ekind
->tcstat
[i
].ekinh
,ekind
->tcstat
[i
].ekinh_old
);
2189 enerd
->term
[F_TEMP
] *= 2; /* result of averages being done over previous and current step,
2190 and there is no previous step */
2192 temp0
= enerd
->term
[F_TEMP
];
2194 /* if using an iterative algorithm, we need to create a working directory for the state. */
2197 bufstate
= init_bufstate(state
);
2201 snew(xcopy
,state
->natoms
);
2202 snew(vcopy
,state
->natoms
);
2203 copy_rvecn(state
->x
,xcopy
,0,state
->natoms
);
2204 copy_rvecn(state
->v
,vcopy
,0,state
->natoms
);
2205 copy_mat(state
->box
,boxcopy
);
2208 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
2209 temperature control */
2210 trotter_seq
= init_npt_vars(ir
,state
,&MassQ
,bTrotter
);
2214 if (constr
&& !ir
->bContinuation
&& ir
->eConstrAlg
== econtLINCS
)
2217 "RMS relative constraint deviation after constraining: %.2e\n",
2218 constr_rmsd(constr
,FALSE
));
2220 fprintf(fplog
,"Initial temperature: %g K\n",enerd
->term
[F_TEMP
]);
2223 fprintf(stderr
,"starting md rerun '%s', reading coordinates from"
2224 " input trajectory '%s'\n\n",
2225 *(top_global
->name
),opt2fn("-rerun",nfile
,fnm
));
2228 fprintf(stderr
,"Calculated time to finish depends on nsteps from "
2229 "run input file,\nwhich may not correspond to the time "
2230 "needed to process input trajectory.\n\n");
2236 fprintf(stderr
,"starting mdrun '%s'\n",
2237 *(top_global
->name
));
2238 if (ir
->nsteps
>= 0)
2240 sprintf(tbuf
,"%8.1f",(ir
->init_step
+ir
->nsteps
)*ir
->delta_t
);
2244 sprintf(tbuf
,"%s","infinite");
2246 if (ir
->init_step
> 0)
2248 fprintf(stderr
,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
2249 gmx_step_str(ir
->init_step
+ir
->nsteps
,sbuf
),tbuf
,
2250 gmx_step_str(ir
->init_step
,sbuf2
),
2251 ir
->init_step
*ir
->delta_t
);
2255 fprintf(stderr
,"%s steps, %s ps.\n",
2256 gmx_step_str(ir
->nsteps
,sbuf
),tbuf
);
2259 fprintf(fplog
,"\n");
2262 /* Set and write start time */
2263 runtime_start(runtime
);
2264 print_date_and_time(fplog
,cr
->nodeid
,"Started mdrun",runtime
);
2265 wallcycle_start(wcycle
,ewcRUN
);
2267 fprintf(fplog
,"\n");
2269 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
2270 /*#ifdef GMX_FAHCORE
2271 chkpt_ret=fcCheckPointParallel( cr->nodeid,
2273 if ( chkpt_ret == 0 )
2274 gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
2278 /***********************************************************
2280 * Loop over MD steps
2282 ************************************************************/
2284 /* if rerunMD then read coordinates and velocities from input trajectory */
2287 if (getenv("GMX_FORCE_UPDATE"))
2289 bForceUpdate
= TRUE
;
2292 bNotLastFrame
= read_first_frame(oenv
,&status
,
2293 opt2fn("-rerun",nfile
,fnm
),
2294 &rerun_fr
,TRX_NEED_X
| TRX_READ_V
);
2295 if (rerun_fr
.natoms
!= top_global
->natoms
)
2298 "Number of atoms in trajectory (%d) does not match the "
2299 "run input file (%d)\n",
2300 rerun_fr
.natoms
,top_global
->natoms
);
2302 if (ir
->ePBC
!= epbcNONE
)
2306 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
);
2308 if (max_cutoff2(ir
->ePBC
,rerun_fr
.box
) < sqr(fr
->rlistlong
))
2310 gmx_fatal(FARGS
,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr
.step
,rerun_fr
.time
);
2313 /* Set the shift vectors.
2314 * Necessary here when have a static box different from the tpr box.
2316 calc_shifts(rerun_fr
.box
,fr
->shift_vec
);
2320 /* loop over MD steps or if rerunMD to end of input trajectory */
2322 /* Skip the first Nose-Hoover integration when we get the state from tpx */
2323 bStateFromTPX
= !opt2bSet("-cpi",nfile
,fnm
);
2324 bInitStep
= bFirstStep
&& (bStateFromTPX
|| bVV
);
2325 bStartingFromCpt
= (Flags
& MD_STARTFROMCPT
) && bInitStep
;
2327 bSumEkinhOld
= FALSE
;
2330 init_global_signals(&gs
,cr
,ir
,repl_ex_nst
);
2332 step
= ir
->init_step
;
2335 if (ir
->nstlist
== -1)
2337 init_nlistheuristics(&nlh
,bGStatEveryStep
,step
);
2340 bLastStep
= (bRerunMD
|| (ir
->nsteps
>= 0 && step_rel
> ir
->nsteps
));
2341 while (!bLastStep
|| (bRerunMD
&& bNotLastFrame
)) {
2343 wallcycle_start(wcycle
,ewcSTEP
);
2345 GMX_MPE_LOG(ev_timestep1
);
2348 if (rerun_fr
.bStep
) {
2349 step
= rerun_fr
.step
;
2350 step_rel
= step
- ir
->init_step
;
2352 if (rerun_fr
.bTime
) {
2362 bLastStep
= (step_rel
== ir
->nsteps
);
2363 t
= t0
+ step
*ir
->delta_t
;
2366 if (ir
->efep
!= efepNO
)
2368 if (bRerunMD
&& rerun_fr
.bLambda
&& (ir
->delta_lambda
!=0))
2370 state_global
->lambda
= rerun_fr
.lambda
;
2374 state_global
->lambda
= lam0
+ step
*ir
->delta_lambda
;
2376 state
->lambda
= state_global
->lambda
;
2377 bDoDHDL
= do_per_step(step
,ir
->nstdhdl
);
2382 update_annealing_target_temp(&(ir
->opts
),t
);
2387 if (!(DOMAINDECOMP(cr
) && !MASTER(cr
)))
2389 for(i
=0; i
<state_global
->natoms
; i
++)
2391 copy_rvec(rerun_fr
.x
[i
],state_global
->x
[i
]);
2395 for(i
=0; i
<state_global
->natoms
; i
++)
2397 copy_rvec(rerun_fr
.v
[i
],state_global
->v
[i
]);
2402 for(i
=0; i
<state_global
->natoms
; i
++)
2404 clear_rvec(state_global
->v
[i
]);
2408 fprintf(stderr
,"\nWARNING: Some frames do not contain velocities.\n"
2409 " Ekin, temperature and pressure are incorrect,\n"
2410 " the virial will be incorrect when constraints are present.\n"
2412 bRerunWarnNoV
= FALSE
;
2416 copy_mat(rerun_fr
.box
,state_global
->box
);
2417 copy_mat(state_global
->box
,state
->box
);
2419 if (vsite
&& (Flags
& MD_RERUN_VSITE
))
2421 if (DOMAINDECOMP(cr
))
2423 gmx_fatal(FARGS
,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
2427 /* Following is necessary because the graph may get out of sync
2428 * with the coordinates if we only have every N'th coordinate set
2430 mk_mshift(fplog
,graph
,fr
->ePBC
,state
->box
,state
->x
);
2431 shift_self(graph
,state
->box
,state
->x
);
2433 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,state
->v
,
2434 top
->idef
.iparams
,top
->idef
.il
,
2435 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
2438 unshift_self(graph
,state
->box
,state
->x
);
2443 /* Stop Center of Mass motion */
2444 bStopCM
= (ir
->comm_mode
!= ecmNO
&& do_per_step(step
,ir
->nstcomm
));
2446 /* Copy back starting coordinates in case we're doing a forcefield scan */
2449 for(ii
=0; (ii
<state
->natoms
); ii
++)
2451 copy_rvec(xcopy
[ii
],state
->x
[ii
]);
2452 copy_rvec(vcopy
[ii
],state
->v
[ii
]);
2454 copy_mat(boxcopy
,state
->box
);
2459 /* for rerun MD always do Neighbour Searching */
2460 bNS
= (bFirstStep
|| ir
->nstlist
!= 0);
2465 /* Determine whether or not to do Neighbour Searching and LR */
2466 bNStList
= (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0);
2468 bNS
= (bFirstStep
|| bExchanged
|| bNStList
||
2469 (ir
->nstlist
== -1 && nlh
.nabnsb
> 0));
2471 if (bNS
&& ir
->nstlist
== -1)
2473 set_nlistheuristics(&nlh
,bFirstStep
|| bExchanged
,step
);
2477 /* < 0 means stop at next step, > 0 means stop at next NS step */
2478 if ( (gs
.set
[eglsSTOPCOND
] < 0 ) ||
2479 ( (gs
.set
[eglsSTOPCOND
] > 0 ) && ( bNS
|| ir
->nstlist
==0)) )
2484 /* Determine whether or not to update the Born radii if doing GB */
2485 bBornRadii
=bFirstStep
;
2486 if (ir
->implicit_solvent
&& (step
% ir
->nstgbradii
==0))
2491 do_log
= do_per_step(step
,ir
->nstlog
) || bFirstStep
|| bLastStep
;
2492 do_verbose
= bVerbose
&&
2493 (step
% stepout
== 0 || bFirstStep
|| bLastStep
);
2495 if (bNS
&& !(bFirstStep
&& ir
->bContinuation
&& !bRerunMD
))
2499 bMasterState
= TRUE
;
2503 bMasterState
= FALSE
;
2504 /* Correct the new box if it is too skewed */
2505 if (DYNAMIC_BOX(*ir
))
2507 if (correct_box(fplog
,step
,state
->box
,graph
))
2509 bMasterState
= TRUE
;
2512 if (DOMAINDECOMP(cr
) && bMasterState
)
2514 dd_collect_state(cr
->dd
,state
,state_global
);
2518 if (DOMAINDECOMP(cr
))
2520 /* Repartition the domain decomposition */
2521 wallcycle_start(wcycle
,ewcDOMDEC
);
2522 dd_partition_system(fplog
,step
,cr
,
2523 bMasterState
,nstglobalcomm
,
2524 state_global
,top_global
,ir
,
2525 state
,&f
,mdatoms
,top
,fr
,
2526 vsite
,shellfc
,constr
,
2527 nrnb
,wcycle
,do_verbose
);
2528 wallcycle_stop(wcycle
,ewcDOMDEC
);
2529 /* If using an iterative integrator, reallocate space to match the decomposition */
2533 if (MASTER(cr
) && do_log
&& !bFFscan
)
2535 print_ebin_header(fplog
,step
,t
,state
->lambda
);
2538 if (ir
->efep
!= efepNO
)
2540 update_mdatoms(mdatoms
,state
->lambda
);
2543 if (bRerunMD
&& rerun_fr
.bV
)
2546 /* We need the kinetic energy at minus the half step for determining
2547 * the full step kinetic energy and possibly for T-coupling.*/
2548 /* This may not be quite working correctly yet . . . . */
2549 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2550 wcycle
,enerd
,NULL
,NULL
,NULL
,NULL
,mu_tot
,
2551 constr
,NULL
,FALSE
,state
->box
,
2552 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2553 CGLO_RERUNMD
| CGLO_GSTAT
| CGLO_TEMPERATURE
);
2555 clear_mat(force_vir
);
2557 /* Ionize the atoms if necessary */
2560 ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
2561 mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
2564 /* Update force field in ffscan program */
2567 if (update_forcefield(fplog,
2569 mdatoms->nr,state->x,state->box)) {
2570 if (gmx_parallel_env_initialized())
2578 GMX_MPE_LOG(ev_timestep2
);
2580 /* We write a checkpoint at this MD step when:
2581 * either at an NS step when we signalled through gs,
2582 * or at the last step (but not when we do not want confout),
2583 * but never at the first step or with rerun.
2585 /* bCPT = (((gs.set[eglsCHKPT] && bNS) ||
2586 (bLastStep && (Flags & MD_CONFOUT))) &&
2587 step > ir->init_step && !bRerunMD);
2590 gs.set[eglsCHKPT] = 0;
2593 /* Determine the energy and pressure:
2594 * at nstcalcenergy steps and at energy output steps (set below).
2596 bNstEner
= (bGStatEveryStep
|| do_per_step(step
,ir
->nstcalcenergy
));
2597 bCalcEnerPres
= bNstEner
;
2599 /* Do we need global communication ? */
2600 bGStat
= (bCalcEnerPres
|| bStopCM
||
2601 (ir
->nstlist
== -1 && !bRerunMD
&& step
>= nlh
.step_nscheck
));
2603 do_ene
= (do_per_step(step
,ir
->nstenergy
) || bLastStep
);
2605 if (do_ene
|| do_log
)
2607 bCalcEnerPres
= TRUE
;
2611 /* these CGLO_ options remain the same throughout the iteration */
2612 cglo_flags
= ((bRerunMD
? CGLO_RERUNMD
: 0) |
2613 (bStopCM
? CGLO_STOPCM
: 0) |
2614 (bGStat
? CGLO_GSTAT
: 0)
2617 force_flags
= (GMX_FORCE_STATECHANGED
|
2618 ((DYNAMIC_BOX(*ir
) || bRerunMD
) ? GMX_FORCE_DYNAMICBOX
: 0) |
2619 GMX_FORCE_ALLFORCES
|
2620 (bNStList
? GMX_FORCE_DOLR
: 0) |
2622 (bCalcEnerPres
? GMX_FORCE_VIRIAL
: 0) |
2623 (bDoDHDL
? GMX_FORCE_DHDL
: 0)
2628 /* Now is the time to relax the shells */
2629 count
=relax_shell_flexcon(fplog
,cr
,bVerbose
,bFFscan
? step
+1 : step
,
2631 bStopCM
,top
,top_global
,
2633 state
,f
,force_vir
,mdatoms
,
2634 nrnb
,wcycle
,graph
,groups
,
2635 shellfc
,fr
,bBornRadii
,t
,mu_tot
,
2636 state
->natoms
,&bConverged
,vsite
,
2647 /* The coordinates (x) are shifted (to get whole molecules)
2649 * This is parallellized as well, and does communication too.
2650 * Check comments in sim_util.c
2653 do_force(fplog
,cr
,ir
,step
,nrnb
,wcycle
,top
,top_global
,groups
,
2654 state
->box
,state
->x
,&state
->hist
,
2655 f
,force_vir
,mdatoms
,enerd
,fcd
,
2656 state
->lambda
,graph
,
2657 fr
,vsite
,mu_tot
,t
,outf
->fp_field
,ed
,bBornRadii
,
2658 (bNS
? GMX_FORCE_NS
: 0) | force_flags
);
2661 GMX_BARRIER(cr
->mpi_comm_mygroup
);
2665 mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
2666 mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
2669 if (bTCR && bFirstStep)
2671 tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
2672 fprintf(fplog,"Done init_coupling\n");
2676 /* ############### START FIRST UPDATE HALF-STEP ############### */
2678 if (bVV
&& !bStartingFromCpt
&& !bRerunMD
)
2684 /* if using velocity verlet with full time step Ekin,
2685 * take the first half step only to compute the
2686 * virial for the first step. From there,
2687 * revert back to the initial coordinates
2688 * so that the input is actually the initial step.
2690 copy_rvecn(state
->v
,cbuf
,0,state
->natoms
); /* should make this better for parallelizing? */
2693 /* this is for NHC in the Ekin(t+dt/2) version of vv */
2696 trotter_update(ir
,step
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
,ettTSEQ2
);
2699 if (ir
->eI
== eiVVAK
)
2701 update_tcouple(fplog
,step
,ir
,state
,ekind
,wcycle
,upd
,&MassQ
,mdatoms
);
2704 update_coords(fplog
,step
,ir
,mdatoms
,state
,
2705 f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
2706 ekind
,M
,wcycle
,upd
,bInitStep
,etrtVELOCITY1
,
2707 cr
,nrnb
,constr
,&top
->idef
);
2711 gmx_iterate_init(&iterate
,bIterations
&& !bInitStep
);
2713 /* for iterations, we save these vectors, as we will be self-consistently iterating
2715 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
2717 /* save the state */
2718 if (bIterations
&& iterate
.bIterate
) {
2719 copy_coupling_state(state
,bufstate
,ekind
,ekind_save
,&(ir
->opts
));
2723 bFirstIterate
= TRUE
;
2724 while (bFirstIterate
|| (bIterations
&& iterate
.bIterate
))
2726 if (bIterations
&& iterate
.bIterate
)
2728 copy_coupling_state(bufstate
,state
,ekind_save
,ekind
,&(ir
->opts
));
2729 if (bFirstIterate
&& bTrotter
)
2731 /* The first time through, we need a decent first estimate
2732 of veta(t+dt) to compute the constraints. Do
2733 this by computing the box volume part of the
2734 trotter integration at this time. Nothing else
2735 should be changed by this routine here. If
2736 !(first time), we start with the previous value
2739 veta_save
= state
->veta
;
2740 trotter_update(ir
,step
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
,trotter_seq
,ettTSEQ0
);
2741 vetanew
= state
->veta
;
2742 state
->veta
= veta_save
;
2747 if ( !bRerunMD
|| rerun_fr
.bV
|| bForceUpdate
) { /* Why is rerun_fr.bV here? Unclear. */
2750 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
2751 &top
->idef
,shake_vir
,NULL
,
2752 cr
,nrnb
,wcycle
,upd
,constr
,
2753 bInitStep
,TRUE
,bCalcEnerPres
,vetanew
);
2755 if (!bOK
&& !bFFscan
)
2757 gmx_fatal(FARGS
,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2762 { /* Need to unshift here if a do_force has been
2763 called in the previous step */
2764 unshift_self(graph
,state
->box
,state
->x
);
2769 /* if VV, compute the pressure and constraints */
2770 /* if VV2, the pressure and constraints only if using pressure control.*/
2771 bPres
= (ir
->eI
==eiVV
|| IR_NPT_TROTTER(ir
));
2772 bTemp
= ((ir
->eI
==eiVV
&&(!bInitStep
)) || (ir
->eI
==eiVVAK
&& IR_NPT_TROTTER(ir
)));
2773 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
2774 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
2775 constr
,NULL
,FALSE
,state
->box
,
2776 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
2779 | (bTemp
? CGLO_TEMPERATURE
:0)
2780 | (bPres
? CGLO_PRESSURE
: 0)
2781 | (bPres
? CGLO_CONSTRAINT
: 0)
2782 | (iterate
.bIterate
? CGLO_ITERATE
: 0)
2783 | (bFirstIterate
? CGLO_FIRSTITERATE
: 0)
2787 /* explanation of above:
2788 a) We compute Ekin at the full time step
2789 if 1) we are using the AveVel Ekin, and it's not the
2790 initial step, or 2) if we are using AveEkin, but need the full
2791 time step kinetic energy for the pressure.
2792 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
2793 EkinAveVel because it's needed for the pressure */
2795 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
2796 if (bVV
&& !bInitStep
)
2798 trotter_update(ir
,step
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
, trotter_seq
,ettTSEQ2
);
2802 done_iterating(cr
,fplog
,step
,&iterate
,bFirstIterate
,
2803 state
->veta
,&vetanew
))
2807 bFirstIterate
= FALSE
;
2810 if (bTrotter
&& !bInitStep
) {
2811 copy_mat(shake_vir
,state
->svir_prev
);
2812 copy_mat(force_vir
,state
->fvir_prev
);
2813 if (IR_NVT_TROTTER(ir
) && ir
->eI
==eiVV
) {
2814 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
2815 enerd
->term
[F_TEMP
] = sum_ekin(&(ir
->opts
),ekind
,NULL
,(ir
->eI
==eiVV
),FALSE
,FALSE
);
2816 enerd
->term
[F_EKIN
] = trace(ekind
->ekin
);
2819 /* if it's the initial step, we performed this first step just to get the constraint virial */
2820 if (bInitStep
&& ir
->eI
==eiVV
) {
2821 copy_rvecn(cbuf
,state
->v
,0,state
->natoms
);
2824 if (fr
->bSepDVDL
&& fplog
&& do_log
)
2826 fprintf(fplog
,sepdvdlformat
,"Constraint",0.0,dvdl
);
2828 enerd
->term
[F_DHDL_CON
] += dvdl
;
2830 GMX_MPE_LOG(ev_timestep1
);
2834 /* MRS -- now done iterating -- compute the conserved quantity */
2837 if (IR_NVT_TROTTER(ir
) || IR_NPT_TROTTER(ir
))
2840 NPT_energy(ir
,state
,&MassQ
);
2841 if ((ir
->eDispCorr
!= edispcEnerPres
) && (ir
->eDispCorr
!= edispcAllEnerPres
))
2843 last_conserved
-= enerd
->term
[F_DISPCORR
];
2847 last_ekin
= enerd
->term
[F_EKIN
]; /* does this get preserved through checkpointing? */
2851 /* ######## END FIRST UPDATE STEP ############## */
2852 /* ######## If doing VV, we now have v(dt) ###### */
2854 /* ################## START TRAJECTORY OUTPUT ################# */
2856 /* Now we have the energies and forces corresponding to the
2857 * coordinates at time t. We must output all of this before
2859 * for RerunMD t is read from input trajectory
2861 GMX_MPE_LOG(ev_output_start
);
2864 if (do_per_step(step
,ir
->nstxout
)) { mdof_flags
|= MDOF_X
; }
2865 if (do_per_step(step
,ir
->nstvout
)) { mdof_flags
|= MDOF_V
; }
2866 if (do_per_step(step
,ir
->nstfout
)) { mdof_flags
|= MDOF_F
; }
2867 if (do_per_step(step
,ir
->nstxtcout
)) { mdof_flags
|= MDOF_XTC
; }
2868 /* if (bCPT) { mdof_flags |= MDOF_CPT; };*/
2872 fcReportProgress( ir
->nsteps
, step
);
2876 /* Enforce writing positions and velocities at end of run */
2877 mdof_flags
|= (MDOF_X
| MDOF_V
);
2879 /* sync bCPT and fc record-keeping */
2880 /* if (bCPT && MASTER(cr))
2881 fcRequestCheckPoint();*/
2884 if (mdof_flags
!= 0)
2886 wallcycle_start(wcycle
,ewcTRAJ
);
2889 if (state->flags & (1<<estLD_RNG))
2891 get_stochd_state(upd,state);
2897 state_global->ekinstate.bUpToDate = FALSE;
2901 update_ekinstate(&state_global->ekinstate,ekind);
2902 state_global->ekinstate.bUpToDate = TRUE;
2904 update_energyhistory(&state_global->enerhist,mdebin);
2907 write_traj(fplog
,cr
,outf
,mdof_flags
,top_global
,
2908 step
,t
,state
,state_global
,f
,f_global
,&n_xtc
,&x_xtc
);
2915 if (bLastStep
&& step_rel
== ir
->nsteps
&&
2916 (Flags
& MD_CONFOUT
) && MASTER(cr
) &&
2917 !bRerunMD
&& !bFFscan
)
2919 /* x and v have been collected in write_traj,
2920 * because a checkpoint file will always be written
2923 fprintf(stderr
,"\nWriting final coordinates.\n");
2924 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
&&
2927 /* Make molecules whole only for confout writing */
2928 do_pbc_mtop(fplog
,ir
->ePBC
,state
->box
,top_global
,state_global
->x
);
2930 /* write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
2931 *top_global->name,top_global,
2932 state_global->x,state_global->v,
2933 ir->ePBC,state->box);*/
2936 wallcycle_stop(wcycle
,ewcTRAJ
);
2938 GMX_MPE_LOG(ev_output_finish
);
2940 /* kludge -- virial is lost with restart for NPT control. Must restart */
2941 if (bStartingFromCpt
&& bVV
)
2943 copy_mat(state
->svir_prev
,shake_vir
);
2944 copy_mat(state
->fvir_prev
,force_vir
);
2946 /* ################## END TRAJECTORY OUTPUT ################ */
2948 /* Determine the pressure:
2949 * always when we want exact averages in the energy file,
2950 * at ns steps when we have pressure coupling,
2951 * otherwise only at energy output steps (set below).
2954 bNstEner
= (bGStatEveryStep
|| do_per_step(step
,ir
->nstcalcenergy
));
2955 bCalcEnerPres
= bNstEner
;
2957 /* Do we need global communication ? */
2958 bGStat
= (bGStatEveryStep
|| bStopCM
|| bNS
||
2959 (ir
->nstlist
== -1 && !bRerunMD
&& step
>= nlh
.step_nscheck
));
2961 do_ene
= (do_per_step(step
,ir
->nstenergy
) || bLastStep
);
2963 if (do_ene
|| do_log
)
2965 bCalcEnerPres
= TRUE
;
2969 /* Determine the wallclock run time up till now */
2970 run_time
= gmx_gettime() - (double)runtime
->real
;
2972 /* Check whether everything is still allright */
2973 if (((int)gmx_get_stop_condition() > handled_stop_condition
)
2979 /* this is just make gs.sig compatible with the hack
2980 of sending signals around by MPI_Reduce with together with
2982 if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns
)
2983 gs
.sig
[eglsSTOPCOND
]=1;
2984 if ( gmx_get_stop_condition() == gmx_stop_cond_next
)
2985 gs
.sig
[eglsSTOPCOND
]=-1;
2986 /* < 0 means stop at next step, > 0 means stop at next NS step */
2990 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2991 gmx_get_signal_name(),
2992 gs
.sig
[eglsSTOPCOND
]==1 ? "NS " : "");
2996 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2997 gmx_get_signal_name(),
2998 gs
.sig
[eglsSTOPCOND
]==1 ? "NS " : "");
3000 handled_stop_condition
=(int)gmx_get_stop_condition();
3002 else if (MASTER(cr
) && (bNS
|| ir
->nstlist
<= 0) &&
3003 (max_hours
> 0 && run_time
> max_hours
*60.0*60.0*0.99) &&
3004 gs
.sig
[eglsSTOPCOND
] == 0 && gs
.set
[eglsSTOPCOND
] == 0)
3006 /* Signal to terminate the run */
3007 gs
.sig
[eglsSTOPCOND
] = 1;
3010 fprintf(fplog
,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step
,sbuf
),max_hours
*0.99);
3012 fprintf(stderr
, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step
,sbuf
),max_hours
*0.99);
3015 if (bResetCountersHalfMaxH
&& MASTER(cr
) &&
3016 run_time
> max_hours
*60.0*60.0*0.495)
3018 gs
.sig
[eglsRESETCOUNTERS
] = 1;
3021 if (ir
->nstlist
== -1 && !bRerunMD
)
3023 /* When bGStatEveryStep=FALSE, global_stat is only called
3024 * when we check the atom displacements, not at NS steps.
3025 * This means that also the bonded interaction count check is not
3026 * performed immediately after NS. Therefore a few MD steps could
3027 * be performed with missing interactions.
3028 * But wrong energies are never written to file,
3029 * since energies are only written after global_stat
3032 if (step
>= nlh
.step_nscheck
)
3034 nlh
.nabnsb
= natoms_beyond_ns_buffer(ir
,fr
,&top
->cgs
,
3035 nlh
.scale_tot
,state
->x
);
3039 /* This is not necessarily true,
3040 * but step_nscheck is determined quite conservatively.
3046 /* In parallel we only have to check for checkpointing in steps
3047 * where we do global communication,
3048 * otherwise the other nodes don't know.
3050 if (MASTER(cr
) && ((bGStat
|| !PAR(cr
)) &&
3053 run_time
>= nchkpt
*cpt_period
*60.0)) &&
3054 gs
.set
[eglsCHKPT
] == 0)
3056 gs
.sig
[eglsCHKPT
] = 1;
3061 gmx_iterate_init(&iterate
,bIterations
);
3064 /* for iterations, we save these vectors, as we will be redoing the calculations */
3065 if (bIterations
&& iterate
.bIterate
)
3067 copy_coupling_state(state
,bufstate
,ekind
,ekind_save
,&(ir
->opts
));
3069 bFirstIterate
= TRUE
;
3070 while (bFirstIterate
|| (bIterations
&& iterate
.bIterate
))
3072 /* We now restore these vectors to redo the calculation with improved extended variables */
3075 copy_coupling_state(bufstate
,state
,ekind_save
,ekind
,&(ir
->opts
));
3078 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
3079 so scroll down for that logic */
3081 /* ######### START SECOND UPDATE STEP ################# */
3082 GMX_MPE_LOG(ev_update_start
);
3084 if (!bRerunMD
|| rerun_fr
.bV
|| bForceUpdate
)
3086 wallcycle_start(wcycle
,ewcUPDATE
);
3088 /* Box is changed in update() when we do pressure coupling,
3089 * but we should still use the old box for energy corrections and when
3090 * writing it to the energy file, so it matches the trajectory files for
3091 * the same timestep above. Make a copy in a separate array.
3093 copy_mat(state
->box
,lastbox
);
3094 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
3097 if (bIterations
&& iterate
.bIterate
)
3105 /* we use a new value of scalevir to converge the iterations faster */
3106 scalevir
= tracevir
/trace(shake_vir
);
3108 msmul(shake_vir
,scalevir
,shake_vir
);
3109 m_add(force_vir
,shake_vir
,total_vir
);
3110 clear_mat(shake_vir
);
3112 trotter_update(ir
,step
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
, trotter_seq
,ettTSEQ3
);
3114 /* We can only do Berendsen coupling after we have summed
3115 * the kinetic energy or virial. Since the happens
3116 * in global_state after update, we should only do it at
3117 * step % nstlist = 1 with bGStatEveryStep=FALSE.
3120 if (ir
->eI
!= eiVVAK
)
3122 update_tcouple(fplog
,step
,ir
,state
,ekind
,wcycle
,upd
,&MassQ
,mdatoms
);
3124 update_pcouple(fplog
,step
,ir
,state
,pcoupl_mu
,M
,wcycle
,
3129 /* velocity half-step update */
3130 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
3131 ekind
,M
,wcycle
,upd
,FALSE
,etrtVELOCITY2
,cr
,nrnb
,constr
,&top
->idef
);
3134 /* Above, initialize just copies ekinh into ekin,
3135 * it doesn't copy position (for VV),
3136 * and entire integrator for MD.
3141 copy_rvecn(state
->x
,cbuf
,0,state
->natoms
);
3144 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
3145 ekind
,M
,wcycle
,upd
,bInitStep
,etrtPOSITION
,cr
,nrnb
,constr
,&top
->idef
);
3146 wallcycle_stop(wcycle
,ewcUPDATE
);
3148 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
3149 &top
->idef
,shake_vir
,force_vir
,
3150 cr
,nrnb
,wcycle
,upd
,constr
,
3151 bInitStep
,FALSE
,bCalcEnerPres
,state
->veta
);
3155 /* erase F_EKIN and F_TEMP here? */
3156 /* just compute the kinetic energy at the half step to perform a trotter step */
3157 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
3158 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
3159 constr
,NULL
,FALSE
,lastbox
,
3160 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
3161 cglo_flags
| CGLO_TEMPERATURE
| CGLO_CONSTRAINT
3163 wallcycle_start(wcycle
,ewcUPDATE
);
3164 trotter_update(ir
,step
,ekind
,enerd
,state
,total_vir
,mdatoms
,&MassQ
, trotter_seq
,ettTSEQ4
);
3165 /* now we know the scaling, we can compute the positions again again */
3166 copy_rvecn(cbuf
,state
->x
,0,state
->natoms
);
3168 update_coords(fplog
,step
,ir
,mdatoms
,state
,f
,fr
->bTwinRange
&& bNStList
,fr
->f_twin
,fcd
,
3169 ekind
,M
,wcycle
,upd
,bInitStep
,etrtPOSITION
,cr
,nrnb
,constr
,&top
->idef
);
3170 wallcycle_stop(wcycle
,ewcUPDATE
);
3172 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
3173 /* are the small terms in the shake_vir here due
3174 * to numerical errors, or are they important
3175 * physically? I'm thinking they are just errors, but not completely sure.
3176 * For now, will call without actually constraining, constr=NULL*/
3177 update_constraints(fplog
,step
,&dvdl
,ir
,ekind
,mdatoms
,state
,graph
,f
,
3178 &top
->idef
,tmp_vir
,force_vir
,
3179 cr
,nrnb
,wcycle
,upd
,NULL
,
3180 bInitStep
,FALSE
,bCalcEnerPres
,state
->veta
);
3182 if (!bOK
&& !bFFscan
)
3184 gmx_fatal(FARGS
,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
3187 if (fr
->bSepDVDL
&& fplog
&& do_log
)
3189 fprintf(fplog
,sepdvdlformat
,"Constraint",0.0,dvdl
);
3191 enerd
->term
[F_DHDL_CON
] += dvdl
;
3195 /* Need to unshift here */
3196 unshift_self(graph
,state
->box
,state
->x
);
3199 GMX_BARRIER(cr
->mpi_comm_mygroup
);
3200 GMX_MPE_LOG(ev_update_finish
);
3204 wallcycle_start(wcycle
,ewcVSITECONSTR
);
3207 shift_self(graph
,state
->box
,state
->x
);
3209 construct_vsites(fplog
,vsite
,state
->x
,nrnb
,ir
->delta_t
,state
->v
,
3210 top
->idef
.iparams
,top
->idef
.il
,
3211 fr
->ePBC
,fr
->bMolPBC
,graph
,cr
,state
->box
);
3215 unshift_self(graph
,state
->box
,state
->x
);
3217 wallcycle_stop(wcycle
,ewcVSITECONSTR
);
3220 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
3221 if (ir
->nstlist
== -1 && bFirstIterate
)
3223 gs
.sig
[eglsNABNSB
] = nlh
.nabnsb
;
3225 compute_globals(fplog
,gstat
,cr
,ir
,fr
,ekind
,state
,state_global
,mdatoms
,nrnb
,vcm
,
3226 wcycle
,enerd
,force_vir
,shake_vir
,total_vir
,pres
,mu_tot
,
3228 bFirstIterate
? &gs
: NULL
,(step
% gs
.nstms
== 0),
3230 top_global
,&pcurr
,top_global
->natoms
,&bSumEkinhOld
,
3232 | (!EI_VV(ir
->eI
) ? CGLO_ENERGY
: 0)
3233 | (!EI_VV(ir
->eI
) ? CGLO_TEMPERATURE
: 0)
3234 | (!EI_VV(ir
->eI
) || bRerunMD
? CGLO_PRESSURE
: 0)
3235 | (bIterations
&& iterate
.bIterate
? CGLO_ITERATE
: 0)
3236 | (bFirstIterate
? CGLO_FIRSTITERATE
: 0)
3239 if (ir
->nstlist
== -1 && bFirstIterate
)
3241 nlh
.nabnsb
= gs
.set
[eglsNABNSB
];
3242 gs
.set
[eglsNABNSB
] = 0;
3244 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
3245 /* ############# END CALC EKIN AND PRESSURE ################# */
3247 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
3248 the virial that should probably be addressed eventually. state->veta has better properies,
3249 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
3250 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
3253 done_iterating(cr
,fplog
,step
,&iterate
,bFirstIterate
,
3254 trace(shake_vir
),&tracevir
))
3258 bFirstIterate
= FALSE
;
3261 update_box(fplog
,step
,ir
,mdatoms
,state
,graph
,f
,
3262 ir
->nstlist
==-1 ? &nlh
.scale_tot
: NULL
,pcoupl_mu
,nrnb
,wcycle
,upd
,bInitStep
,FALSE
);
3264 /* ################# END UPDATE STEP 2 ################# */
3265 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
3267 /* The coordinates (x) were unshifted in update */
3268 /* if (bFFscan && (shellfc==NULL || bConverged))
3270 if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
3272 &(top_global->mols),mdatoms->massT,pres))
3274 if (gmx_parallel_env_initialized())
3278 fprintf(stderr,"\n");
3284 /* We will not sum ekinh_old,
3285 * so signal that we still have to do it.
3287 bSumEkinhOld
= TRUE
;
3292 /* Only do GCT when the relaxation of shells (minimization) has converged,
3293 * otherwise we might be coupling to bogus energies.
3294 * In parallel we must always do this, because the other sims might
3298 /* Since this is called with the new coordinates state->x, I assume
3299 * we want the new box state->box too. / EL 20040121
3301 /* do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
3303 mdatoms,&(top->idef),mu_aver,
3304 top_global->mols.nr,cr,
3305 state->box,total_vir,pres,
3306 mu_tot,state->x,f,bConverged);
3310 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
3312 sum_dhdl(enerd
,state
->lambda
,ir
);
3313 /* use the directly determined last velocity, not actually the averaged half steps */
3314 if (bTrotter
&& ir
->eI
==eiVV
)
3316 enerd
->term
[F_EKIN
] = last_ekin
;
3318 enerd
->term
[F_ETOT
] = enerd
->term
[F_EPOT
] + enerd
->term
[F_EKIN
];
3327 if (IR_NVT_TROTTER(ir
)) {
3328 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] + last_conserved
;
3330 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] +
3331 NPT_energy(ir
,state
,&MassQ
);
3335 enerd
->term
[F_ECONSERVED
] =
3336 enerd
->term
[F_ETOT
] + vrescale_energy(&(ir
->opts
),
3337 state
->therm_integral
);
3343 /* Check for excessively large energies */
3347 real etot_max = 1e200;
3349 real etot_max = 1e30;
3351 if (fabs(enerd->term[F_ETOT]) > etot_max)
3353 fprintf(stderr,"Energy too large (%g), giving up\n",
3354 enerd->term[F_ETOT]);
3357 /* ######### END PREPARING EDR OUTPUT ########### */
3359 /* Time for performance */
3360 if (((step
% stepout
) == 0) || bLastStep
)
3362 runtime_upd_proc(runtime
);
3368 gmx_bool do_dr
,do_or
;
3370 if (!(bStartingFromCpt
&& (EI_VV(ir
->eI
))))
3374 upd_mdebin(mdebin
,bDoDHDL
,TRUE
,
3375 t
,mdatoms
->tmass
,enerd
,state
,lastbox
,
3376 shake_vir
,force_vir
,total_vir
,pres
,
3377 ekind
,mu_tot
,constr
);
3381 upd_mdebin_step(mdebin
);
3384 do_dr
= do_per_step(step
,ir
->nstdisreout
);
3385 do_or
= do_per_step(step
,ir
->nstorireout
);
3387 print_ebin(outf
->fp_ene
,do_ene
,do_dr
,do_or
,do_log
?fplog
:NULL
,
3389 eprNORMAL
,bCompact
,mdebin
,fcd
,groups
,&(ir
->opts
));
3391 if (ir
->ePull
!= epullNO
)
3393 pull_print_output(ir
->pull
,step
,t
);
3396 if (do_per_step(step
,ir
->nstlog
))
3398 if(fflush(fplog
) != 0)
3400 gmx_fatal(FARGS
,"Cannot flush logfile - maybe you are out of quota?");
3406 /* Remaining runtime */
3407 if (MULTIMASTER(cr
) && (do_verbose
|| gmx_got_usr_signal() ))
3411 fprintf(stderr
,"\n");
3413 print_time(stderr
,runtime
,step
,ir
,cr
);
3416 /* Set new positions for the group to embed */
3422 } else if (step_rel
<=(it_xy
+it_z
))
3426 resize(ins_at
,r_ins
,state_global
->x
,pos_ins
,fac
);
3429 /* Replica exchange */
3430 /* bExchanged = FALSE;
3431 if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
3432 do_per_step(step,repl_ex_nst))
3434 bExchanged = replica_exchange(fplog,cr,repl_ex,
3435 state_global,enerd->term,
3438 if (bExchanged && PAR(cr))
3440 if (DOMAINDECOMP(cr))
3442 dd_partition_system(fplog,step,cr,TRUE,1,
3443 state_global,top_global,ir,
3444 state,&f,mdatoms,top,fr,
3445 vsite,shellfc,constr,
3450 bcast_state(cr,state,FALSE);
3456 bStartingFromCpt
= FALSE
;
3458 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
3459 /* Complicated conditional when bGStatEveryStep=FALSE.
3460 * We can not just use bGStat, since then the simulation results
3461 * would depend on nstenergy and nstlog or step_nscheck.
3463 if (((state
->flags
& (1<<estPRES_PREV
)) ||
3464 (state
->flags
& (1<<estSVIR_PREV
)) ||
3465 (state
->flags
& (1<<estFVIR_PREV
))) &&
3467 (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0) ||
3468 (ir
->nstlist
< 0 && nlh
.nabnsb
> 0) ||
3469 (ir
->nstlist
== 0 && bGStat
)))
3471 /* Store the pressure in t_state for pressure coupling
3472 * at the next MD step.
3474 if (state
->flags
& (1<<estPRES_PREV
))
3476 copy_mat(pres
,state
->pres_prev
);
3480 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
3484 /* read next frame from input trajectory */
3485 bNotLastFrame
= read_next_frame(oenv
,status
,&rerun_fr
);
3488 if (!bRerunMD
|| !rerun_fr
.bStep
)
3490 /* increase the MD step number */
3495 cycles
= wallcycle_stop(wcycle
,ewcSTEP
);
3496 if (DOMAINDECOMP(cr
) && wcycle
)
3498 dd_cycles_add(cr
->dd
,cycles
,ddCyclStep
);
3501 if (step_rel
== wcycle_get_reset_counters(wcycle
) ||
3502 gs
.set
[eglsRESETCOUNTERS
] != 0)
3504 /* Reset all the counters related to performance over the run */
3505 reset_all_counters(fplog
,cr
,step
,&step_rel
,ir
,wcycle
,nrnb
,runtime
);
3506 wcycle_set_reset_counters(wcycle
,-1);
3507 bResetCountersHalfMaxH
= FALSE
;
3508 gs
.set
[eglsRESETCOUNTERS
] = 0;
3511 /* End of main MD loop */
3513 write_sto_conf_mtop(ftp2fn(efSTO
,nfile
,fnm
),
3514 *top_global
->name
,top_global
,
3515 state_global
->x
,state_global
->v
,
3516 ir
->ePBC
,state
->box
);
3519 runtime_end(runtime
);
3526 if (!(cr
->duty
& DUTY_PME
))
3528 /* Tell the PME only node to finish */
3534 if (ir
->nstcalcenergy
> 0 && !bRerunMD
)
3536 print_ebin(outf
->fp_ene
,FALSE
,FALSE
,FALSE
,fplog
,step
,t
,
3537 eprAVER
,FALSE
,mdebin
,fcd
,groups
,&(ir
->opts
));
3545 if (ir
->nstlist
== -1 && nlh
.nns
> 0 && fplog
)
3547 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
)));
3548 fprintf(fplog
,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh
.ab
/nlh
.nns
);
3551 if (shellfc
&& fplog
)
3553 fprintf(fplog
,"Fraction of iterations that converged: %.2f %%\n",
3554 (nconverged
*100.0)/step_rel
);
3555 fprintf(fplog
,"Average number of force evaluations per MD step: %.2f\n\n",
3559 /* if (repl_ex_nst > 0 && MASTER(cr))
3561 print_replica_exchange_statistics(fplog,repl_ex);
3564 runtime
->nsteps_done
= step_rel
;
3570 int mdrunner_membed(FILE *fplog
,t_commrec
*cr
,int nfile
,const t_filenm fnm
[],
3571 const output_env_t oenv
, gmx_bool bVerbose
,gmx_bool bCompact
,
3573 ivec ddxyz
,int dd_node_order
,real rdd
,real rconstr
,
3574 const char *dddlb_opt
,real dlb_scale
,
3575 const char *ddcsx
,const char *ddcsy
,const char *ddcsz
,
3576 int nstepout
,int resetstep
,int nmultisim
,int repl_ex_nst
,int repl_ex_seed
,
3577 real pforce
,real cpt_period
,real max_hours
,
3578 const char *deviceOptions
,
3579 unsigned long Flags
,
3580 real xy_fac
, real xy_max
, real z_fac
, real z_max
,
3581 int it_xy
, int it_z
, real probe_rad
, int low_up_rm
,
3582 int pieces
, gmx_bool bALLOW_ASYMMETRY
, int maxwarn
)
3584 double nodetime
=0,realtime
;
3585 t_inputrec
*inputrec
;
3586 t_state
*state
=NULL
;
3589 int npme_major
,npme_minor
;
3592 gmx_mtop_t
*mtop
=NULL
;
3593 t_mdatoms
*mdatoms
=NULL
;
3594 t_forcerec
*fr
=NULL
;
3597 gmx_pme_t
*pmedata
=NULL
;
3598 gmx_vsite_t
*vsite
=NULL
;
3599 gmx_constr_t constr
;
3600 int i
,m
,nChargePerturbed
=-1,status
,nalloc
;
3602 gmx_wallcycle_t wcycle
;
3603 gmx_bool bReadRNG
,bReadEkin
;
3605 gmx_runtime_t runtime
;
3607 gmx_large_int_t reset_counters
;
3608 gmx_edsam_t ed
=NULL
;
3609 t_commrec
*cr_old
=cr
;
3610 int nthreads
=1,nthreads_requested
=1;
3614 int rm_bonded_at
,fr_id
,fr_i
=0,tmp_id
,warn
=0;
3615 int ng
,j
,max_lip_rm
,ins_grp_id
,ins_nat
,mem_nat
,ntype
,lip_rm
,tpr_version
;
3616 real xy_step
=0,z_step
=0;
3618 rvec
*r_ins
=NULL
,fac
;
3619 t_block
*ins_at
,*rest_at
;
3623 gmx_groups_t
*groups
;
3624 gmx_bool bExcl
=FALSE
;
3627 char **piecename
=NULL
;
3629 /* CAUTION: threads may be started later on in this function, so
3630 cr doesn't reflect the final parallel state right now */
3634 if (bVerbose
&& SIMMASTER(cr
))
3636 fprintf(stderr
,"Getting Loaded...\n");
3639 if (Flags
& MD_APPENDFILES
)
3647 /* Read (nearly) all data required for the simulation */
3648 read_tpx_state(ftp2fn(efTPX
,nfile
,fnm
),inputrec
,state
,NULL
,mtop
);
3650 /* NOW the threads will be started: */
3654 /* END OF CAUTION: cr is now reliable */
3658 /* now broadcast everything to the non-master nodes/threads: */
3659 init_parallel(fplog
, cr
, inputrec
, mtop
);
3661 /* now make sure the state is initialized and propagated */
3662 set_state_entries(state
,inputrec
,cr
->nnodes
);
3664 if (can_use_allvsall(inputrec
,mtop
,TRUE
,cr
,fplog
))
3666 /* All-vs-all loops do not work with domain decomposition */
3667 Flags
|= MD_PARTDEC
;
3670 if (!EEL_PME(inputrec
->coulombtype
) || (Flags
& MD_PARTDEC
))
3679 tpr_version
= get_tpr_version(ftp2fn(efTPX
,nfile
,fnm
));
3681 gmx_fatal(FARGS
,"Version of *.tpr file to old (%d). Rerun grompp with gromacs VERSION 4.0.3 or newer.\n",tpr_version
);
3683 if( inputrec
->eI
!= eiMD
)
3684 gmx_input("Change integrator to md in mdp file.");
3687 gmx_input("Sorry, parallel g_membed is not yet fully functrional.");
3689 groups
=&(mtop
->groups
);
3691 atoms
=gmx_mtop_global_atoms(mtop
);
3693 fprintf(stderr
,"\nSelect a group to embed in the membrane:\n");
3694 get_index(&atoms
,ftp2fn_null(efNDX
,nfile
,fnm
),1,&(ins_at
->nr
),&(ins_at
->index
),&ins
);
3695 ins_grp_id
= search_string(ins
,groups
->ngrpname
,(groups
->grpname
));
3696 fprintf(stderr
,"\nSelect a group to embed %s into (e.g. the membrane):\n",ins
);
3697 get_index(&atoms
,ftp2fn_null(efNDX
,nfile
,fnm
),1,&(mem_p
->mem_at
.nr
),&(mem_p
->mem_at
.index
),&(mem_p
->name
));
3699 pos_ins
->pieces
=pieces
;
3700 snew(pos_ins
->nidx
,pieces
);
3701 snew(pos_ins
->subindex
,pieces
);
3702 snew(piecename
,pieces
);
3705 fprintf(stderr
,"\nSelect pieces to embed:\n");
3706 get_index(&atoms
,ftp2fn_null(efNDX
,nfile
,fnm
),pieces
,pos_ins
->nidx
,pos_ins
->subindex
,piecename
);
3710 /*use whole embedded group*/
3711 snew(pos_ins
->nidx
,1);
3712 snew(pos_ins
->subindex
,1);
3713 pos_ins
->nidx
[0]=ins_at
->nr
;
3714 pos_ins
->subindex
[0]=ins_at
->index
;
3717 if(probe_rad
<0.2199999)
3720 fprintf(stderr
,"\nWarning %d:\nA probe radius (-rad) smaller than 0.2 can result in overlap between waters "
3721 "and the group to embed, which will result in Lincs errors etc.\nIf you are sure, you can increase maxwarn.\n\n",warn
);
3724 if(xy_fac
<0.09999999)
3727 fprintf(stderr
,"\nWarning %d:\nThe initial size of %s is probably too smal.\n"
3728 "If you are sure, you can increase maxwarn.\n\n",warn
,ins
);
3734 fprintf(stderr
,"\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d) is probably too small.\n"
3735 "Increase -nxy or, if you are sure, you can increase maxwarn.\n\n",warn
,ins
,it_xy
);
3738 if( (it_z
<100) && ( z_fac
<0.99999999 || z_fac
>1.0000001) )
3741 fprintf(stderr
,"\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d) is probably too small.\n"
3742 "Increase -nz or, if you are sure, you can increase maxwarn.\n\n",warn
,ins
,it_z
);
3745 if(it_xy
+it_z
>inputrec
->nsteps
)
3748 fprintf(stderr
,"\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the number of steps in the tpr.\n"
3749 "If you are sure, you can increase maxwarn.\n\n",warn
);
3753 if( inputrec
->opts
.ngfrz
==1)
3754 gmx_fatal(FARGS
,"You did not specify \"%s\" as a freezegroup.",ins
);
3755 for(i
=0;i
<inputrec
->opts
.ngfrz
;i
++)
3757 tmp_id
= mtop
->groups
.grps
[egcFREEZE
].nm_ind
[i
];
3758 if(ins_grp_id
==tmp_id
)
3765 gmx_fatal(FARGS
,"\"%s\" not as freezegroup defined in the mdp-file.",ins
);
3768 if( inputrec
->opts
.nFreeze
[fr_i
][i
] != 1)
3769 gmx_fatal(FARGS
,"freeze dimensions for %s are not Y Y Y\n",ins
);
3771 ng
= groups
->grps
[egcENER
].nr
;
3773 gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
3779 if (inputrec
->opts
.egp_flags
[ng
*i
+j
] == EGP_EXCL
)
3782 if ( (groups
->grps
[egcENER
].nm_ind
[i
] != ins_grp_id
) || (groups
->grps
[egcENER
].nm_ind
[j
] != ins_grp_id
) )
3783 gmx_fatal(FARGS
,"Energy exclusions \"%s\" and \"%s\" do not match the group to embed \"%s\"",
3784 *groups
->grpname
[groups
->grps
[egcENER
].nm_ind
[i
]],
3785 *groups
->grpname
[groups
->grps
[egcENER
].nm_ind
[j
]],ins
);
3790 gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in the freeze group");
3792 /* Set all atoms in box*/
3793 /*set_inbox(state->natoms,state->x);*/
3795 /* Guess the area the protein will occupy in the membrane plane Calculate area per lipid*/
3797 ins_nat
= init_ins_at(ins_at
,rest_at
,state
,pos_ins
,groups
,ins_grp_id
,xy_max
);
3798 /* Check moleculetypes in insertion group */
3799 check_types(ins_at
,rest_at
,mtop
);
3801 mem_nat
= init_mem_at(mem_p
,mtop
,state
->x
,state
->box
,pos_ins
);
3803 prot_area
= est_prot_area(pos_ins
,state
->x
,ins_at
,mem_p
);
3804 if ( (prot_area
>7.5) && ( (state
->box
[XX
][XX
]*state
->box
[YY
][YY
]-state
->box
[XX
][YY
]*state
->box
[YY
][XX
])<50) )
3807 fprintf(stderr
,"\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
3808 "This might cause pressure problems during the growth phase. Just try with\n"
3809 "current setup (-maxwarn + 1), but if pressure problems occur, lower the\n"
3810 "compressibility in the mdp-file or use no pressure coupling at all.\n\n",warn
);
3813 gmx_fatal(FARGS
,"Too many warnings.\n");
3815 printf("The estimated area of the protein in the membrane is %.3f nm^2\n",prot_area
);
3816 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
);
3818 /* Maximum number of lipids to be removed*/
3819 max_lip_rm
=(int)(2*prot_area
/mem_p
->lip_area
);
3820 printf("Maximum number of lipids that will be removed is %d.\n",max_lip_rm
);
3822 printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
3823 "This resizing will be done with respect to the geometrical center of all protein atoms\n"
3824 "that span the membrane region, i.e. z between %.3f and %.3f\n\n",xy_fac
,z_fac
,mem_p
->zmin
,mem_p
->zmax
);
3826 /* resize the protein by xy and by z if necessary*/
3827 snew(r_ins
,ins_at
->nr
);
3828 init_resize(ins_at
,r_ins
,pos_ins
,mem_p
,state
->x
,bALLOW_ASYMMETRY
);
3829 fac
[0]=fac
[1]=xy_fac
;
3832 xy_step
=(xy_max
-xy_fac
)/(double)(it_xy
);
3833 z_step
=(z_max
-z_fac
)/(double)(it_z
-1);
3835 resize(ins_at
,r_ins
,state
->x
,pos_ins
,fac
);
3837 /* remove overlapping lipids and water from the membrane box*/
3838 /*mark molecules to be removed*/
3840 set_pbc(pbc
,inputrec
->ePBC
,state
->box
);
3843 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
);
3844 lip_rm
-= low_up_rm
;
3847 for(i
=0;i
<rm_p
->nr
;i
++)
3848 fprintf(fplog
,"rm mol %d\n",rm_p
->mol
[i
]);
3850 for(i
=0;i
<mtop
->nmolblock
;i
++)
3853 for(j
=0;j
<rm_p
->nr
;j
++)
3854 if(rm_p
->block
[j
]==i
)
3856 printf("Will remove %d %s molecules\n",ntype
,*(mtop
->moltype
[mtop
->molblock
[i
].type
].name
));
3859 if(lip_rm
>max_lip_rm
)
3862 fprintf(stderr
,"\nWarning %d:\nTrying to remove a larger lipid area than the estimated protein area\n"
3863 "Try making the -xyinit resize factor smaller. If you are sure about this increase maxwarn.\n\n",warn
);
3866 /*remove all lipids and waters overlapping and update all important structures*/
3867 rm_group(inputrec
,groups
,mtop
,rm_p
,state
,ins_at
,pos_ins
);
3869 rm_bonded_at
= rm_bonded(ins_at
,mtop
);
3870 if (rm_bonded_at
!= ins_at
->nr
)
3872 fprintf(stderr
,"Warning: The number of atoms for which the bonded interactions are removed is %d, "
3873 "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
3874 "molecule type as atoms that are not to be embedded.\n",rm_bonded_at
,ins_at
->nr
);
3878 gmx_fatal(FARGS
,"Too many warnings.\nIf you are sure these warnings are harmless, you can increase -maxwarn");
3882 if (ftp2bSet(efTOP
,nfile
,fnm
))
3883 top_update(opt2fn("-p",nfile
,fnm
),ins
,rm_p
,mtop
);
3891 fcRegisterSteps(inputrec
->nsteps
,inputrec
->init_step
);
3894 /* NMR restraints must be initialized before load_checkpoint,
3895 * since with time averaging the history is added to t_state.
3896 * For proper consistency check we therefore need to extend
3898 * So the PME-only nodes (if present) will also initialize
3899 * the distance restraints.
3903 /* This needs to be called before read_checkpoint to extend the state */
3904 init_disres(fplog
,mtop
,inputrec
,cr
,Flags
& MD_PARTDEC
,fcd
,state
);
3906 if (gmx_mtop_ftype_count(mtop
,F_ORIRES
) > 0)
3908 if (PAR(cr
) && !(Flags
& MD_PARTDEC
))
3910 gmx_fatal(FARGS
,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
3912 /* Orientation restraints */
3915 init_orires(fplog
,mtop
,state
->x
,inputrec
,cr
->ms
,&(fcd
->orires
),
3920 if (DEFORM(*inputrec
))
3922 /* Store the deform reference box before reading the checkpoint */
3925 copy_mat(state
->box
,box
);
3929 gmx_bcast(sizeof(box
),box
,cr
);
3931 /* Because we do not have the update struct available yet
3932 * in which the reference values should be stored,
3933 * we store them temporarily in static variables.
3934 * This should be thread safe, since they are only written once
3935 * and with identical values.
3937 /* deform_init_init_step_tpx = inputrec->init_step;*/
3938 /* copy_mat(box,deform_init_box_tpx);*/
3941 if (opt2bSet("-cpi",nfile
,fnm
))
3943 /* Check if checkpoint file exists before doing continuation.
3944 * This way we can use identical input options for the first and subsequent runs...
3946 if( gmx_fexist_master(opt2fn_master("-cpi",nfile
,fnm
,cr
),cr
) )
3948 load_checkpoint(opt2fn_master("-cpi",nfile
,fnm
,cr
),&fplog
,
3949 cr
,Flags
& MD_PARTDEC
,ddxyz
,
3950 inputrec
,state
,&bReadRNG
,&bReadEkin
,
3951 (Flags
& MD_APPENDFILES
));
3955 Flags
|= MD_READ_RNG
;
3959 Flags
|= MD_READ_EKIN
;
3964 if ((MASTER(cr
) || (Flags
& MD_SEPPOT
)) && (Flags
& MD_APPENDFILES
))
3966 gmx_log_open(ftp2fn(efLOG
,nfile
,fnm
),cr
,!(Flags
& MD_SEPPOT
),
3972 copy_mat(state
->box
,box
);
3977 gmx_bcast(sizeof(box
),box
,cr
);
3980 if (bVerbose
&& SIMMASTER(cr
))
3982 fprintf(stderr
,"Loaded with Money\n\n");
3985 if (PAR(cr
) && !((Flags
& MD_PARTDEC
) || EI_TPI(inputrec
->eI
)))
3987 cr
->dd
= init_domain_decomposition(fplog
,cr
,Flags
,ddxyz
,rdd
,rconstr
,
3988 dddlb_opt
,dlb_scale
,
3992 &ddbox
,&npme_major
,&npme_minor
);
3994 make_dd_communicators(fplog
,cr
,dd_node_order
);
3996 /* Set overallocation to avoid frequent reallocation of arrays */
3997 set_over_alloc_dd(TRUE
);
4001 /* PME, if used, is done on all nodes with 1D decomposition */
4003 cr
->duty
= (DUTY_PP
| DUTY_PME
);
4004 npme_major
= cr
->nnodes
;
4007 if (inputrec
->ePBC
== epbcSCREW
)
4010 "pbc=%s is only implemented with domain decomposition",
4011 epbc_names
[inputrec
->ePBC
]);
4017 /* After possible communicator splitting in make_dd_communicators.
4018 * we can set up the intra/inter node communication.
4020 gmx_setup_nodecomm(fplog
,cr
);
4023 wcycle
= wallcycle_init(fplog
,resetstep
,cr
);
4026 /* Master synchronizes its value of reset_counters with all nodes
4027 * including PME only nodes */
4028 reset_counters
= wcycle_get_reset_counters(wcycle
);
4029 gmx_bcast_sim(sizeof(reset_counters
),&reset_counters
,cr
);
4030 wcycle_set_reset_counters(wcycle
, reset_counters
);
4035 if (cr
->duty
& DUTY_PP
)
4037 /* For domain decomposition we allocate dynamically
4038 * in dd_partition_system.
4040 if (DOMAINDECOMP(cr
))
4042 bcast_state_setup(cr
,state
);
4052 bcast_state(cr
,state
,TRUE
);
4056 /* Dihedral Restraints */
4057 if (gmx_mtop_ftype_count(mtop
,F_DIHRES
) > 0)
4059 init_dihres(fplog
,mtop
,inputrec
,fcd
);
4062 /* Initiate forcerecord */
4064 init_forcerec(fplog
,oenv
,fr
,fcd
,inputrec
,mtop
,cr
,box
,FALSE
,
4065 opt2fn("-table",nfile
,fnm
),
4066 opt2fn("-tablep",nfile
,fnm
),
4067 opt2fn("-tableb",nfile
,fnm
),FALSE
,pforce
);
4069 /* version for PCA_NOT_READ_NODE (see md.c) */
4070 /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
4071 "nofile","nofile","nofile",FALSE,pforce);
4073 fr
->bSepDVDL
= ((Flags
& MD_SEPPOT
) == MD_SEPPOT
);
4075 /* Initialize QM-MM */
4078 init_QMMMrec(cr
,box
,mtop
,inputrec
,fr
);
4081 /* Initialize the mdatoms structure.
4082 * mdatoms is not filled with atom data,
4083 * as this can not be done now with domain decomposition.
4085 mdatoms
= init_mdatoms(fplog
,mtop
,inputrec
->efep
!=efepNO
);
4087 /* Initialize the virtual site communication */
4088 vsite
= init_vsite(mtop
,cr
);
4090 calc_shifts(box
,fr
->shift_vec
);
4092 /* With periodic molecules the charge groups should be whole at start up
4093 * and the virtual sites should not be far from their proper positions.
4095 if (!inputrec
->bContinuation
&& MASTER(cr
) &&
4096 !(inputrec
->ePBC
!= epbcNONE
&& inputrec
->bPeriodicMols
))
4098 /* Make molecules whole at start of run */
4099 if (fr
->ePBC
!= epbcNONE
)
4101 do_pbc_first_mtop(fplog
,inputrec
->ePBC
,box
,mtop
,state
->x
);
4105 /* Correct initial vsite positions are required
4106 * for the initial distribution in the domain decomposition
4107 * and for the initial shell prediction.
4109 construct_vsites_mtop(fplog
,vsite
,mtop
,state
->x
);
4113 /* Initiate PPPM if necessary */
4114 if (fr
->eeltype
== eelPPPM
)
4116 if (mdatoms
->nChargePerturbed
)
4118 gmx_fatal(FARGS
,"Free energy with %s is not implemented",
4119 eel_names
[fr
->eeltype
]);
4121 status
= gmx_pppm_init(fplog
,cr
,oenv
,FALSE
,TRUE
,box
,
4122 getenv("GMXGHAT"),inputrec
, (Flags
& MD_REPRODUCIBLE
));
4125 gmx_fatal(FARGS
,"Error %d initializing PPPM",status
);
4129 if (EEL_PME(fr
->eeltype
))
4131 ewaldcoeff
= fr
->ewaldcoeff
;
4132 pmedata
= &fr
->pmedata
;
4141 /* This is a PME only node */
4143 /* We don't need the state */
4146 ewaldcoeff
= calc_ewaldcoeff(inputrec
->rcoulomb
, inputrec
->ewald_rtol
);
4150 /* Initiate PME if necessary,
4151 * either on all nodes or on dedicated PME nodes only. */
4152 if (EEL_PME(inputrec
->coulombtype
))
4156 nChargePerturbed
= mdatoms
->nChargePerturbed
;
4158 if (cr
->npmenodes
> 0)
4160 /* The PME only nodes need to know nChargePerturbed */
4161 gmx_bcast_sim(sizeof(nChargePerturbed
),&nChargePerturbed
,cr
);
4163 if (cr
->duty
& DUTY_PME
)
4165 status
= gmx_pme_init(pmedata
,cr
,npme_major
,npme_minor
,inputrec
,
4166 mtop
? mtop
->natoms
: 0,nChargePerturbed
,
4167 (Flags
& MD_REPRODUCIBLE
));
4170 gmx_fatal(FARGS
,"Error %d initializing PME",status
);
4176 /* if (integrator[inputrec->eI].func == do_md
4179 integrator[inputrec->eI].func == do_md_openmm
4183 /* Turn on signal handling on all nodes */
4185 * (A user signal from the PME nodes (if any)
4186 * is communicated to the PP nodes.
4188 signal_handler_install();
4191 if (cr
->duty
& DUTY_PP
)
4193 if (inputrec
->ePull
!= epullNO
)
4195 /* Initialize pull code */
4196 init_pull(fplog
,inputrec
,nfile
,fnm
,mtop
,cr
,oenv
,
4197 EI_DYNAMICS(inputrec
->eI
) && MASTER(cr
),Flags
);
4200 constr
= init_constraints(fplog
,mtop
,inputrec
,ed
,state
,cr
);
4202 if (DOMAINDECOMP(cr
))
4204 dd_init_bondeds(fplog
,cr
->dd
,mtop
,vsite
,constr
,inputrec
,
4205 Flags
& MD_DDBONDCHECK
,fr
->cginfo_mb
);
4207 set_dd_parameters(fplog
,cr
->dd
,dlb_scale
,inputrec
,fr
,&ddbox
);
4209 setup_dd_grid(fplog
,cr
->dd
);
4212 /* Now do whatever the user wants us to do (how flexible...) */
4213 do_md_membed(fplog
,cr
,nfile
,fnm
,
4214 oenv
,bVerbose
,bCompact
,
4217 nstepout
,inputrec
,mtop
,
4219 mdatoms
,nrnb
,wcycle
,ed
,fr
,
4220 repl_ex_nst
,repl_ex_seed
,
4221 cpt_period
,max_hours
,
4225 fac
, r_ins
, pos_ins
, ins_at
,
4226 xy_step
, z_step
, it_xy
, it_z
);
4228 if (inputrec
->ePull
!= epullNO
)
4230 finish_pull(fplog
,inputrec
->pull
);
4236 gmx_pmeonly(*pmedata
,cr
,nrnb
,wcycle
,ewaldcoeff
,FALSE
,inputrec
);
4239 if (EI_DYNAMICS(inputrec
->eI
) || EI_TPI(inputrec
->eI
))
4241 /* Some timing stats */
4244 if (runtime
.proc
== 0)
4246 runtime
.proc
= runtime
.real
;
4255 wallcycle_stop(wcycle
,ewcRUN
);
4257 /* Finish up, write some stuff
4258 * if rerunMD, don't write last frame again
4260 finish_run(fplog
,cr
,ftp2fn(efSTO
,nfile
,fnm
),
4261 inputrec
,nrnb
,wcycle
,&runtime
,
4262 EI_DYNAMICS(inputrec
->eI
) && !MULTISIM(cr
));
4264 /* Does what it says */
4265 print_date_and_time(fplog
,cr
->nodeid
,"Finished mdrun",&runtime
);
4267 /* Close logfile already here if we were appending to it */
4268 if (MASTER(cr
) && (Flags
& MD_APPENDFILES
))
4270 gmx_log_close(fplog
);
4278 rc
=(int)gmx_get_stop_condition();
4283 int gmx_membed(int argc
,char *argv
[])
4285 const char *desc
[] = {
4286 "g_membed embeds a membrane protein into an equilibrated lipid bilayer at the position",
4287 "and orientation specified by the user.\n",
4289 "SHORT MANUAL\n------------\n",
4290 "The user should merge the structure files of the protein and membrane (+solvent), creating a",
4291 "single structure file with the protein overlapping the membrane at the desired position and",
4292 "orientation. Box size should be taken from the membrane structure file. The corresponding topology",
4293 "files should also be merged. Consecutively, create a tpr file (input for g_membed) from these files,"
4294 "with the following options included in the mdp file.\n",
4295 " - integrator = md\n",
4296 " - energygrp = Protein (or other group that you want to insert)\n",
4297 " - freezegrps = Protein\n",
4298 " - freezedim = Y Y Y\n",
4299 " - energygrp_excl = Protein Protein\n",
4300 "The output is a structure file containing the protein embedded in the membrane. If a topology",
4301 "file is provided, the number of lipid and ",
4302 "solvent molecules will be updated to match the new structure file.\n",
4303 "For a more extensive manual see Wolf et al, J Comp Chem 31 (2010) 2169-2174, Appendix.\n",
4305 "SHORT METHOD DESCRIPTION\n",
4306 "------------------------\n",
4307 "1. The protein is resized around its center of mass by a factor -xy in the xy-plane",
4308 "(the membrane plane) and a factor -z in the z-direction (if the size of the",
4309 "protein in the z-direction is the same or smaller than the width of the membrane, a",
4310 "-z value larger than 1 can prevent that the protein will be enveloped by the lipids).\n",
4311 "2. All lipid and solvent molecules overlapping with the resized protein are removed. All",
4312 "intraprotein interactions are turned off to prevent numerical issues for small values of -xy",
4314 "3. One md step is performed.\n",
4315 "4. The resize factor (-xy or -z) is incremented by a small amount ((1-xy)/nxy or (1-z)/nz) and the",
4316 "protein is resized again around its center of mass. The resize factor for the xy-plane",
4317 "is incremented first. The resize factor for the z-direction is not changed until the -xy factor",
4318 "is 1 (thus after -nxy iteration).\n",
4319 "5. Repeat step 3 and 4 until the protein reaches its original size (-nxy + -nz iterations).\n",
4320 "For a more extensive method descrition see Wolf et al, J Comp Chem, 31 (2010) 2169-2174.\n",
4323 " - Protein can be any molecule you want to insert in the membrane.\n",
4324 " - It is recommended to perform a short equilibration run after the embedding",
4325 "(see Wolf et al, J Comp Chem 31 (2010) 2169-2174, to re-equilibrate the membrane. Clearly",
4326 "protein equilibration might require longer.\n",
4331 { efTPX
, "-f", "into_mem", ffREAD
},
4332 { efNDX
, "-n", "index", ffOPTRD
},
4333 { efTOP
, "-p", "topol", ffOPTRW
},
4334 { efTRN
, "-o", NULL
, ffWRITE
},
4335 { efXTC
, "-x", NULL
, ffOPTWR
},
4336 { efCPT
, "-cpi", NULL
, ffOPTRD
},
4337 { efCPT
, "-cpo", NULL
, ffOPTWR
},
4338 { efSTO
, "-c", "membedded", ffWRITE
},
4339 { efEDR
, "-e", "ener", ffWRITE
},
4340 { efLOG
, "-g", "md", ffWRITE
},
4341 { efEDI
, "-ei", "sam", ffOPTRD
},
4342 { efTRX
, "-rerun", "rerun", ffOPTRD
},
4343 { efXVG
, "-table", "table", ffOPTRD
},
4344 { efXVG
, "-tablep", "tablep", ffOPTRD
},
4345 { efXVG
, "-tableb", "table", ffOPTRD
},
4346 { efXVG
, "-dhdl", "dhdl", ffOPTWR
},
4347 { efXVG
, "-field", "field", ffOPTWR
},
4348 { efXVG
, "-table", "table", ffOPTRD
},
4349 { efXVG
, "-tablep", "tablep", ffOPTRD
},
4350 { efXVG
, "-tableb", "table", ffOPTRD
},
4351 { efTRX
, "-rerun", "rerun", ffOPTRD
},
4352 { efXVG
, "-tpi", "tpi", ffOPTWR
},
4353 { efXVG
, "-tpid", "tpidist", ffOPTWR
},
4354 { efEDI
, "-ei", "sam", ffOPTRD
},
4355 { efEDO
, "-eo", "sam", ffOPTWR
},
4356 { efGCT
, "-j", "wham", ffOPTRD
},
4357 { efGCT
, "-jo", "bam", ffOPTWR
},
4358 { efXVG
, "-ffout", "gct", ffOPTWR
},
4359 { efXVG
, "-devout", "deviatie", ffOPTWR
},
4360 { efXVG
, "-runav", "runaver", ffOPTWR
},
4361 { efXVG
, "-px", "pullx", ffOPTWR
},
4362 { efXVG
, "-pf", "pullf", ffOPTWR
},
4363 { efMTX
, "-mtx", "nm", ffOPTWR
},
4364 { efNDX
, "-dn", "dipole", ffOPTWR
}
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 -rdd > 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 -rerun" },
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 -c 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 -maxh" },
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
];
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();
4553 if (repl_ex_nst
!= 0 && nmultisim
< 2)
4554 gmx_fatal(FARGS
,"Need at least two replicas for replica exchange (option -multi)");
4556 if (nmultisim
> 1) {
4558 init_multisystem(cr
,nmultisim
,NFILE
,fnm
,TRUE
);
4560 gmx_fatal(FARGS
,"mdrun -multi is not supported with the thread library.Please compile GROMACS with MPI support");
4564 /* Check if there is ANY checkpoint file available */
4566 sim_part_fn
= sim_part
;
4567 if (opt2bSet("-cpi",NFILE
,fnm
))
4570 read_checkpoint_simulation_part(opt2fn_master("-cpi", NFILE
,fnm
,cr
),
4571 &sim_part_fn
,NULL
,cr
,
4572 bAppendFiles
,NFILE
,fnm
,
4573 part_suffix
,&bAddPart
);
4574 if (sim_part_fn
==0 && MASTER(cr
))
4576 fprintf(stdout
,"No previous checkpoint file present, assuming this is a new run.\n");
4580 sim_part
= sim_part_fn
+ 1;
4585 bAppendFiles
= FALSE
;
4590 sim_part_fn
= sim_part
;
4593 if (bAddPart
&& sim_part_fn
> 1)
4595 /* This is a continuation run, rename trajectory output files
4596 (except checkpoint files) */
4597 /* create new part name first (zero-filled) */
4598 sprintf(suffix
,"%s%04d",part_suffix
,sim_part_fn
);
4600 add_suffix_to_output_names(fnm
,NFILE
,suffix
);
4601 fprintf(stdout
,"Checkpoint file is from part %d, new output files will be suffixed '%s'.\n",sim_part
-1,suffix
);
4604 Flags
= opt2bSet("-rerun",NFILE
,fnm
) ? MD_RERUN
: 0;
4605 Flags
= Flags
| (bSepPot
? MD_SEPPOT
: 0);
4606 Flags
= Flags
| (bIonize
? MD_IONIZE
: 0);
4607 Flags
= Flags
| (bPartDec
? MD_PARTDEC
: 0);
4608 Flags
= Flags
| (bDDBondCheck
? MD_DDBONDCHECK
: 0);
4609 Flags
= Flags
| (bDDBondComm
? MD_DDBONDCOMM
: 0);
4610 Flags
= Flags
| (bConfout
? MD_CONFOUT
: 0);
4611 Flags
= Flags
| (bRerunVSite
? MD_RERUN_VSITE
: 0);
4612 Flags
= Flags
| (bReproducible
? MD_REPRODUCIBLE
: 0);
4613 Flags
= Flags
| (bAppendFiles
? MD_APPENDFILES
: 0);
4614 Flags
= Flags
| (sim_part
>1 ? MD_STARTFROMCPT
: 0);
4615 Flags
= Flags
| (bResetCountersHalfWay
? MD_RESETCOUNTERSHALFWAY
: 0);
4618 /* We postpone opening the log file if we are appending, so we can
4619 first truncate the old log file and append to the correct position
4621 if ((MASTER(cr
) || bSepPot
) && !bAppendFiles
)
4623 gmx_log_open(ftp2fn(efLOG
,NFILE
,fnm
),cr
,!bSepPot
,Flags
,&fplog
);
4624 CopyRight(fplog
,argv
[0]);
4625 please_cite(fplog
,"Hess2008b");
4626 please_cite(fplog
,"Spoel2005a");
4627 please_cite(fplog
,"Lindahl2001a");
4628 please_cite(fplog
,"Berendsen95a");
4635 ddxyz
[XX
] = (int)(realddxyz
[XX
] + 0.5);
4636 ddxyz
[YY
] = (int)(realddxyz
[YY
] + 0.5);
4637 ddxyz
[ZZ
] = (int)(realddxyz
[ZZ
] + 0.5);
4639 /* even if nthreads = 1, we still call this one */
4641 rc
= mdrunner_membed(fplog
, cr
, NFILE
, fnm
, oenv
, bVerbose
, bCompact
,
4643 ddxyz
, dd_node_order
, rdd
, rconstr
, dddlb_opt
[0], dlb_scale
,
4644 ddcsx
, ddcsy
, ddcsz
, nstepout
, resetstep
, nmultisim
, repl_ex_nst
,
4645 repl_ex_seed
, pforce
, cpt_period
, max_hours
, deviceOptions
, Flags
,
4646 xy_fac
,xy_max
,z_fac
,z_max
,
4647 it_xy
,it_z
,probe_rad
,low_up_rm
,
4648 pieces
,bALLOW_ASYMMETRY
,maxwarn
);
4650 if (gmx_parallel_env_initialized())
4653 if (MULTIMASTER(cr
)) {
4657 /* Log file has to be closed in mdrunner if we are appending to it
4658 (fplog not set here) */
4659 fprintf(stderr
,"Please cite:\nWolf et al, J Comp Chem 31 (2010) 2169-2174.\n");
4661 if (MASTER(cr
) && !bAppendFiles
)
4663 gmx_log_close(fplog
);