?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Updated code checks with distance restraints
[gromacs.git] / src / tools / gmx_membed.c
blobdb2854f40eafdf52a1e43fcf153ea448d0b5712a
1 /*
2 *
3 * This source code is part of
4 *
5 * G R O M A C S
6 *
7 * GROningen MAchine for Chemical Simulations
8 *
9 * VERSION 3.2.0
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
14
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
19 *
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
26 *
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
29 *
30 * For more info, check our website at http://www.gromacs.org
31 *
32 * And Hey:
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
34 */
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38
39 #include <signal.h>
40 #include <stdlib.h>
41 #include "typedefs.h"
42 #include "smalloc.h"
43 #include "sysstuff.h"
44 #include "vec.h"
45 #include "statutil.h"
46 #include "macros.h"
47 #include "copyrite.h"
48 #include "main.h"
49 #include "futil.h"
50 #include "edsam.h"
51 #include "checkpoint.h"
52 #include "vcm.h"
53 #include "mdebin.h"
54 #include "nrnb.h"
55 #include "calcmu.h"
56 #include "index.h"
57 #include "vsite.h"
58 #include "update.h"
59 #include "ns.h"
60 #include "trnio.h"
61 #include "xtcio.h"
62 #include "mdrun.h"
63 #include "confio.h"
64 #include "network.h"
65 #include "pull.h"
66 #include "xvgr.h"
67 #include "physics.h"
68 #include "names.h"
69 #include "disre.h"
70 #include "orires.h"
71 #include "dihre.h"
72 #include "pppm.h"
73 #include "pme.h"
74 #include "mdatoms.h"
75 #include "qmmm.h"
76 #include "mpelogging.h"
77 #include "domdec.h"
78 #include "partdec.h"
79 #include "topsort.h"
80 #include "coulomb.h"
81 #include "constr.h"
82 #include "shellfc.h"
83 #include "mvdata.h"
84 #include "checkpoint.h"
85 #include "mtop_util.h"
86 #include "tpxio.h"
87 #include "string2.h"
88 #include "sighandler.h"
89 #include "gmx_ana.h"
90
91 #ifdef GMX_LIB_MPI
92 #include <mpi.h>
93 #endif
94 #ifdef GMX_THREADS
95 #include "tmpi.h"
96 #endif
97
98 /* afm stuf */
99 #include "pull.h"
100
101 /* We use the same defines as in mvdata.c here */
102 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d),(cr))
103 #define nblock_bc(cr,nr,d) gmx_bcast((nr)*sizeof((d)[0]), (d),(cr))
104 #define snew_bc(cr,d,nr) { if (!MASTER(cr)) snew((d),(nr)); }
105
106 /* The following two variables and the signal_handler function
107 * is used from pme.c as well
108 */
109
110 typedef struct {
111 t_state s;
112 rvec *f;
113 real epot;
114 real fnorm;
115 real fmax;
116 int a_fmax;
117 } em_state_t;
118
119 typedef struct {
120 int it_xy;
121 int it_z;
122 int xy_step;
123 int z_step;
124 rvec xmin;
125 rvec xmax;
126 rvec *geom_cent;
127 int pieces;
128 int *nidx;
129 atom_id **subindex;
130 } pos_ins_t;
131
132 typedef struct {
133 int id;
134 char *name;
135 int nr;
136 int natoms; /*nr of atoms per lipid*/
137 int mol1; /*id of the first lipid molecule*/
138 real area;
139 } lip_t;
140
141 typedef struct {
142 char *name;
143 t_block mem_at;
144 int *mol_id;
145 int nmol;
146 real lip_area;
147 real zmin;
148 real zmax;
149 real zmed;
150 } mem_t;
151
152 typedef struct {
153 int *mol;
154 int *block;
155 int nr;
156 } rmm_t;
157
158 int search_string(char *s,int ng,char ***gn)
159 {
160 int i;
161
162 for(i=0; (i<ng); i++)
163 if (gmx_strcasecmp(s,*gn[i]) == 0)
164 return i;
165
166 gmx_fatal(FARGS,"Group %s not found in indexfile.\nMaybe you have non-default groups in your mdp file, while not using the '-n' option of grompp.\nIn that case use the '-n' option.\n",s);
167
168 return -1;
169 }
170
171 int get_mol_id(int at,int nmblock,gmx_molblock_t *mblock, int *type, int *block)
172 {
173 int mol_id=0;
174 int i;
175
176 for(i=0;i<nmblock;i++)
177 {
178 if(at<(mblock[i].nmol*mblock[i].natoms_mol))
179 {
180 mol_id+=at/mblock[i].natoms_mol;
181 *type = mblock[i].type;
182 *block = i;
183 return mol_id;
184 } else {
185 at-= mblock[i].nmol*mblock[i].natoms_mol;
186 mol_id+=mblock[i].nmol;
187 }
188 }
189
190 gmx_fatal(FARGS,"Something is wrong in mol ids, at %d, mol_id %d",at,mol_id);
191
192 return -1;
193 }
194
195 int get_block(int mol_id,int nmblock,gmx_molblock_t *mblock)
196 {
197 int i;
198 int nmol=0;
199
200 for(i=0;i<nmblock;i++)
201 {
202 nmol+=mblock[i].nmol;
203 if(mol_id<nmol)
204 return i;
205 }
206
207 gmx_fatal(FARGS,"mol_id %d larger than total number of molecules %d.\n",mol_id,nmol);
208
209 return -1;
210 }
211
212 int get_tpr_version(const char *infile)
213 {
214 char buf[STRLEN];
215 gmx_bool bDouble;
216 int precision,fver;
217 t_fileio *fio;
218
219 fio = open_tpx(infile,"r");
220 gmx_fio_checktype(fio);
221
222 precision = sizeof(real);
223
224 gmx_fio_do_string(fio,buf);
225 if (strncmp(buf,"VERSION",7))
226 gmx_fatal(FARGS,"Can not read file %s,\n"
227 " this file is from a Gromacs version which is older than 2.0\n"
228 " Make a new one with grompp or use a gro or pdb file, if possible",
229 gmx_fio_getname(fio));
230 gmx_fio_do_int(fio,precision);
231 bDouble = (precision == sizeof(double));
232 if ((precision != sizeof(float)) && !bDouble)
233 gmx_fatal(FARGS,"Unknown precision in file %s: real is %d bytes "
234 "instead of %d or %d",
235 gmx_fio_getname(fio),precision,sizeof(float),sizeof(double));
236 gmx_fio_setprecision(fio,bDouble);
237 fprintf(stderr,"Reading file %s, %s (%s precision)\n",
238 gmx_fio_getname(fio),buf,bDouble ? "double" : "single");
239
240 gmx_fio_do_int(fio,fver);
241
242 close_tpx(fio);
243
244 return fver;
245 }
246
247 void set_inbox(int natom, rvec *x)
248 {
249 rvec tmp;
250 int i;
251
252 tmp[XX]=tmp[YY]=tmp[ZZ]=0.0;
253 for(i=0;i<natom;i++)
254 {
255 if(x[i][XX]<tmp[XX]) tmp[XX]=x[i][XX];
256 if(x[i][YY]<tmp[YY]) tmp[YY]=x[i][YY];
257 if(x[i][ZZ]<tmp[ZZ]) tmp[ZZ]=x[i][ZZ];
258 }
259
260 for(i=0;i<natom;i++)
261 rvec_inc(x[i],tmp);
262 }
263
264 int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist)
265 {
266 int i,j,nr,mol_id;
267 int type=0,block=0;
268 gmx_bool bNEW;
269
270 nr=0;
271 snew(tlist->index,at->nr);
272 for (i=0;i<at->nr;i++)
273 {
274 bNEW=TRUE;
275 mol_id = get_mol_id(at->index[i],mtop->nmolblock,mtop->molblock,&type,&block);
276 for(j=0;j<nr;j++)
277 {
278 if(tlist->index[j]==type)
279 bNEW=FALSE;
280 }
281 if(bNEW==TRUE)
282 {
283 tlist->index[nr]=type;
284 nr++;
285 }
286 }
287
288 srenew(tlist->index,nr);
289 return nr;
290 }
291
292 void check_types(t_block *ins_at,t_block *rest_at,gmx_mtop_t *mtop)
293 {
294 t_block *ins_mtype,*rest_mtype;
295 int i,j;
296
297 snew(ins_mtype,1);
298 snew(rest_mtype,1);
299 ins_mtype->nr = get_mtype_list(ins_at , mtop, ins_mtype );
300 rest_mtype->nr = get_mtype_list(rest_at, mtop, rest_mtype);
301
302 for(i=0;i<ins_mtype->nr;i++)
303 {
304 for(j=0;j<rest_mtype->nr;j++)
305 {
306 if(ins_mtype->index[i]==rest_mtype->index[j])
307 gmx_fatal(FARGS,"Moleculetype %s is found both in the group to insert and the rest of the system.\n"
308 "Because we need to exclude all interactions between the atoms in the group to\n"
309 "insert, the same moleculetype can not be used in both groups. Change the\n"
310 "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n"
311 "an appropriate *.itp file",*(mtop->moltype[rest_mtype->index[j]].name),
312 *(mtop->moltype[rest_mtype->index[j]].name));
313 }
314 }
315
316 sfree(ins_mtype->index);
317 sfree(rest_mtype->index);
318 sfree(ins_mtype);
319 sfree(rest_mtype);
320 }
321
322 int init_ins_at(t_block *ins_at,t_block *rest_at,t_state *state, pos_ins_t *pos_ins,gmx_groups_t *groups,int ins_grp_id, real xy_max)
323 {
324 int i,gid,c=0;
325 real x,xmin,xmax,y,ymin,ymax,z,zmin,zmax;
326
327 snew(rest_at->index,state->natoms);
328
329 xmin=xmax=state->x[ins_at->index[0]][XX];
330 ymin=ymax=state->x[ins_at->index[0]][YY];
331 zmin=zmax=state->x[ins_at->index[0]][ZZ];
332
333 for(i=0;i<state->natoms;i++)
334 {
335 gid = groups->grpnr[egcFREEZE][i];
336 if(groups->grps[egcFREEZE].nm_ind[gid]==ins_grp_id)
337 {
338 x=state->x[i][XX];
339 if (x<xmin) xmin=x;
340 if (x>xmax) xmax=x;
341 y=state->x[i][YY];
342 if (y<ymin) ymin=y;
343 if (y>ymax) ymax=y;
344 z=state->x[i][ZZ];
345 if (z<zmin) zmin=z;
346 if (z>zmax) zmax=z;
347 } else {
348 rest_at->index[c]=i;
349 c++;
350 }
351 }
352
353 rest_at->nr=c;
354 srenew(rest_at->index,c);
355
356 if(xy_max>1.000001)
357 {
358 pos_ins->xmin[XX]=xmin-((xmax-xmin)*xy_max-(xmax-xmin))/2;
359 pos_ins->xmin[YY]=ymin-((ymax-ymin)*xy_max-(ymax-ymin))/2;
360
361 pos_ins->xmax[XX]=xmax+((xmax-xmin)*xy_max-(xmax-xmin))/2;
362 pos_ins->xmax[YY]=ymax+((ymax-ymin)*xy_max-(ymax-ymin))/2;
363 } else {
364 pos_ins->xmin[XX]=xmin;
365 pos_ins->xmin[YY]=ymin;
366
367 pos_ins->xmax[XX]=xmax;
368 pos_ins->xmax[YY]=ymax;
369 }
370
371 /* 6.0 is estimated thickness of bilayer */
372 if( (zmax-zmin) < 6.0 )
373 {
374 pos_ins->xmin[ZZ]=zmin+(zmax-zmin)/2.0-3.0;
375 pos_ins->xmax[ZZ]=zmin+(zmax-zmin)/2.0+3.0;
376 } else {
377 pos_ins->xmin[ZZ]=zmin;
378 pos_ins->xmax[ZZ]=zmax;
379 }
380
381 return c;
382 }
383
384 real est_prot_area(pos_ins_t *pos_ins,rvec *r,t_block *ins_at, mem_t *mem_p)
385 {
386 real x,y,dx=0.15,dy=0.15,area=0.0;
387 real add;
388 int c,at;
389
390 for(x=pos_ins->xmin[XX];x<pos_ins->xmax[XX];x+=dx)
391 {
392 for(y=pos_ins->xmin[YY];y<pos_ins->xmax[YY];y+=dy)
393 {
394 c=0;
395 add=0.0;
396 do
397 {
398 at=ins_at->index[c];
399 if ( (r[at][XX]>=x) && (r[at][XX]<x+dx) &&
400 (r[at][YY]>=y) && (r[at][YY]<y+dy) &&
401 (r[at][ZZ]>mem_p->zmin+1.0) && (r[at][ZZ]<mem_p->zmax-1.0) )
402 add=1.0;
403 c++;
404 } while ( (c<ins_at->nr) && (add<0.5) );
405 area+=add;
406 }
407 }
408 area=area*dx*dy;
409
410 return area;
411 }
412
413 void init_lip(matrix box, gmx_mtop_t *mtop, lip_t *lip)
414 {
415 int i;
416 real mem_area;
417 int mol1=0;
418
419 mem_area = box[XX][XX]*box[YY][YY]-box[XX][YY]*box[YY][XX];
420 for(i=0;i<mtop->nmolblock;i++)
421 {
422 if(mtop->molblock[i].type == lip->id)
423 {
424 lip->nr=mtop->molblock[i].nmol;
425 lip->natoms=mtop->molblock[i].natoms_mol;
426 }
427 }
428 lip->area=2.0*mem_area/(double)lip->nr;
429
430 for (i=0;i<lip->id;i++)
431 mol1+=mtop->molblock[i].nmol;
432 lip->mol1=mol1;
433 }
434
435 int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins)
436 {
437 int i,j,at,mol,nmol,nmolbox,count;
438 t_block *mem_a;
439 real z,zmin,zmax,mem_area;
440 gmx_bool bNew;
441 atom_id *mol_id;
442 int type=0,block=0;
443
444 nmol=count=0;
445 mem_a=&(mem_p->mem_at);
446 snew(mol_id,mem_a->nr);
447 /* snew(index,mem_a->nr); */
448 zmin=pos_ins->xmax[ZZ];
449 zmax=pos_ins->xmin[ZZ];
450 for(i=0;i<mem_a->nr;i++)
451 {
452 at=mem_a->index[i];
453 if( (r[at][XX]>pos_ins->xmin[XX]) && (r[at][XX]<pos_ins->xmax[XX]) &&
454 (r[at][YY]>pos_ins->xmin[YY]) && (r[at][YY]<pos_ins->xmax[YY]) &&
455 (r[at][ZZ]>pos_ins->xmin[ZZ]) && (r[at][ZZ]<pos_ins->xmax[ZZ]) )
456 {
457 mol = get_mol_id(at,mtop->nmolblock,mtop->molblock,&type,&block);
458
459 bNew=TRUE;
460 for(j=0;j<nmol;j++)
461 if(mol == mol_id[j])
462 bNew=FALSE;
463
464 if(bNew)
465 {
466 mol_id[nmol]=mol;
467 nmol++;
468 }
469
470 z=r[at][ZZ];
471 if(z<zmin) zmin=z;
472 if(z>zmax) zmax=z;
473
474 /* index[count]=at;*/
475 count++;
476 }
477 }
478
479 mem_p->nmol=nmol;
480 srenew(mol_id,nmol);
481 mem_p->mol_id=mol_id;
482 /* srenew(index,count);*/
483 /* mem_p->mem_at.nr=count;*/
484 /* sfree(mem_p->mem_at.index);*/
485 /* mem_p->mem_at.index=index;*/
486
487 if((zmax-zmin)>(box[ZZ][ZZ]-0.5))
488 gmx_fatal(FARGS,"Something is wrong with your membrane. Max and min z values are %f and %f.\n"
489 "Maybe your membrane is not centered in the box, but located at the box edge in the z-direction,\n"
490 "so that one membrane is distributed over two periodic box images. Another possibility is that\n"
491 "your water layer is not thick enough.\n",zmax,zmin);
492 mem_p->zmin=zmin;
493 mem_p->zmax=zmax;
494 mem_p->zmed=(zmax-zmin)/2+zmin;
495
496 /*number of membrane molecules in protein box*/
497 nmolbox = count/mtop->molblock[block].natoms_mol;
498 /*mem_area = box[XX][XX]*box[YY][YY]-box[XX][YY]*box[YY][XX];
499 mem_p->lip_area = 2.0*mem_area/(double)mem_p->nmol;*/
500 mem_area = (pos_ins->xmax[XX]-pos_ins->xmin[XX])*(pos_ins->xmax[YY]-pos_ins->xmin[YY]);
501 mem_p->lip_area = 2.0*mem_area/(double)nmolbox;
502
503 return mem_p->mem_at.nr;
504 }
505
506 void init_resize(t_block *ins_at,rvec *r_ins,pos_ins_t *pos_ins,mem_t *mem_p,rvec *r, gmx_bool bALLOW_ASYMMETRY)
507 {
508 int i,j,at,c,outsidesum,gctr=0;
509 int idxsum=0;
510
511 /*sanity check*/
512 for (i=0;i<pos_ins->pieces;i++)
513 idxsum+=pos_ins->nidx[i];
514 if (idxsum!=ins_at->nr)
515 gmx_fatal(FARGS,"Piecewise sum of inserted atoms not same as size of group selected to insert.");
516
517 snew(pos_ins->geom_cent,pos_ins->pieces);
518 for (i=0;i<pos_ins->pieces;i++)
519 {
520 c=0;
521 outsidesum=0;
522 for(j=0;j<DIM;j++)
523 pos_ins->geom_cent[i][j]=0;
524
525 for(j=0;j<DIM;j++)
526 pos_ins->geom_cent[i][j]=0;
527 for (j=0;j<pos_ins->nidx[i];j++)
528 {
529 at=pos_ins->subindex[i][j];
530 copy_rvec(r[at],r_ins[gctr]);
531 if( (r_ins[gctr][ZZ]<mem_p->zmax) && (r_ins[gctr][ZZ]>mem_p->zmin) )
532 {
533 rvec_inc(pos_ins->geom_cent[i],r_ins[gctr]);
534 c++;
535 }
536 else
537 outsidesum++;
538 gctr++;
539 }
540 if (c>0)
541 svmul(1/(double)c,pos_ins->geom_cent[i],pos_ins->geom_cent[i]);
542 if (!bALLOW_ASYMMETRY)
543 pos_ins->geom_cent[i][ZZ]=mem_p->zmed;
544
545 fprintf(stderr,"Embedding piece %d with center of geometry: %f %f %f\n",i,pos_ins->geom_cent[i][XX],pos_ins->geom_cent[i][YY],pos_ins->geom_cent[i][ZZ]);
546 }
547 fprintf(stderr,"\n");
548 }
549
550 void resize(t_block *ins_at, rvec *r_ins, rvec *r, pos_ins_t *pos_ins,rvec fac)
551 {
552 int i,j,k,at,c=0;
553 for (k=0;k<pos_ins->pieces;k++)
554 for(i=0;i<pos_ins->nidx[k];i++)
555 {
556 at=pos_ins->subindex[k][i];
557 for(j=0;j<DIM;j++)
558 r[at][j]=pos_ins->geom_cent[k][j]+fac[j]*(r_ins[c][j]-pos_ins->geom_cent[k][j]);
559 c++;
560 }
561 }
562
563 int gen_rm_list(rmm_t *rm_p,t_block *ins_at,t_block *rest_at,t_pbc *pbc, gmx_mtop_t *mtop,
564 rvec *r, rvec *r_ins, mem_t *mem_p, pos_ins_t *pos_ins, real probe_rad, int low_up_rm, gmx_bool bALLOW_ASYMMETRY)
565 {
566 int i,j,k,l,at,at2,mol_id;
567 int type=0,block=0;
568 int nrm,nupper,nlower;
569 real r_min_rad,z_lip,min_norm;
570 gmx_bool bRM;
571 rvec dr,dr_tmp;
572 real *dist;
573 int *order;
574
575 r_min_rad=probe_rad*probe_rad;
576 snew(rm_p->mol,mtop->mols.nr);
577 snew(rm_p->block,mtop->mols.nr);
578 nrm=nupper=0;
579 nlower=low_up_rm;
580 for(i=0;i<ins_at->nr;i++)
581 {
582 at=ins_at->index[i];
583 for(j=0;j<rest_at->nr;j++)
584 {
585 at2=rest_at->index[j];
586 pbc_dx(pbc,r[at],r[at2],dr);
587
588 if(norm2(dr)<r_min_rad)
589 {
590 mol_id = get_mol_id(at2,mtop->nmolblock,mtop->molblock,&type,&block);
591 bRM=TRUE;
592 for(l=0;l<nrm;l++)
593 if(rm_p->mol[l]==mol_id)
594 bRM=FALSE;
595 if(bRM)
596 {
597 /*fprintf(stderr,"%d wordt toegevoegd\n",mol_id);*/
598 rm_p->mol[nrm]=mol_id;
599 rm_p->block[nrm]=block;
600 nrm++;
601 z_lip=0.0;
602 for(l=0;l<mem_p->nmol;l++)
603 {
604 if(mol_id==mem_p->mol_id[l])
605 {
606 for(k=mtop->mols.index[mol_id];k<mtop->mols.index[mol_id+1];k++)
607 z_lip+=r[k][ZZ];
608 z_lip/=mtop->molblock[block].natoms_mol;
609 if(z_lip<mem_p->zmed)
610 nlower++;
611 else
612 nupper++;
613 }
614 }
615 }
616 }
617 }
618 }
619
620 /*make sure equal number of lipids from upper and lower layer are removed */
621 if( (nupper!=nlower) && (!bALLOW_ASYMMETRY) )
622 {
623 snew(dist,mem_p->nmol);
624 snew(order,mem_p->nmol);
625 for(i=0;i<mem_p->nmol;i++)
626 {
627 at = mtop->mols.index[mem_p->mol_id[i]];
628 pbc_dx(pbc,r[at],pos_ins->geom_cent[0],dr);
629 if (pos_ins->pieces>1)
630 {
631 /*minimum dr value*/
632 min_norm=norm2(dr);
633 for (k=1;k<pos_ins->pieces;k++)
634 {
635 pbc_dx(pbc,r[at],pos_ins->geom_cent[k],dr_tmp);
636 if (norm2(dr_tmp) < min_norm)
637 {
638 min_norm=norm2(dr_tmp);
639 copy_rvec(dr_tmp,dr);
640 }
641 }
642 }
643 dist[i]=dr[XX]*dr[XX]+dr[YY]*dr[YY];
644 j=i-1;
645 while (j>=0 && dist[i]<dist[order[j]])
646 {
647 order[j+1]=order[j];
648 j--;
649 }
650 order[j+1]=i;
651 }
652
653 i=0;
654 while(nupper!=nlower)
655 {
656 mol_id=mem_p->mol_id[order[i]];
657 block=get_block(mol_id,mtop->nmolblock,mtop->molblock);
658
659 bRM=TRUE;
660 for(l=0;l<nrm;l++)
661 if(rm_p->mol[l]==mol_id)
662 bRM=FALSE;
663 if(bRM)
664 {
665 z_lip=0;
666 for(k=mtop->mols.index[mol_id];k<mtop->mols.index[mol_id+1];k++)
667 z_lip+=r[k][ZZ];
668 z_lip/=mtop->molblock[block].natoms_mol;
669 if(nupper>nlower && z_lip<mem_p->zmed)
670 {
671 rm_p->mol[nrm]=mol_id;
672 rm_p->block[nrm]=block;
673 nrm++;
674 nlower++;
675 }
676 else if (nupper<nlower && z_lip>mem_p->zmed)
677 {
678 rm_p->mol[nrm]=mol_id;
679 rm_p->block[nrm]=block;
680 nrm++;
681 nupper++;
682 }
683 }
684 i++;
685
686 if(i>mem_p->nmol)
687 gmx_fatal(FARGS,"Trying to remove more lipid molecules than there are in the membrane");
688 }
689 sfree(dist);
690 sfree(order);
691 }
692
693 rm_p->nr=nrm;
694 srenew(rm_p->mol,nrm);
695 srenew(rm_p->block,nrm);
696
697 return nupper+nlower;
698 }
699
700 void rm_group(t_inputrec *ir, gmx_groups_t *groups, gmx_mtop_t *mtop, rmm_t *rm_p, t_state *state, t_block *ins_at, pos_ins_t *pos_ins)
701 {
702 int i,j,k,n,rm,mol_id,at,block;
703 rvec *x_tmp,*v_tmp;
704 atom_id *list,*new_mols;
705 unsigned char *new_egrp[egcNR];
706 gmx_bool bRM;
707
708 snew(list,state->natoms);
709 n=0;
710 for(i=0;i<rm_p->nr;i++)
711 {
712 mol_id=rm_p->mol[i];
713 at=mtop->mols.index[mol_id];
714 block =rm_p->block[i];
715 mtop->molblock[block].nmol--;
716 for(j=0;j<mtop->molblock[block].natoms_mol;j++)
717 {
718 list[n]=at+j;
719 n++;
720 }
721
722 mtop->mols.index[mol_id]=-1;
723 }
724
725 mtop->mols.nr-=rm_p->nr;
726 mtop->mols.nalloc_index-=rm_p->nr;
727 snew(new_mols,mtop->mols.nr);
728 for(i=0;i<mtop->mols.nr+rm_p->nr;i++)
729 {
730 j=0;
731 if(mtop->mols.index[i]!=-1)
732 {
733 new_mols[j]=mtop->mols.index[i];
734 j++;
735 }
736 }
737 sfree(mtop->mols.index);
738 mtop->mols.index=new_mols;
739
740
741 mtop->natoms-=n;
742 state->natoms-=n;
743 state->nalloc=state->natoms;
744 snew(x_tmp,state->nalloc);
745 snew(v_tmp,state->nalloc);
746
747 for(i=0;i<egcNR;i++)
748 {
749 if(groups->grpnr[i]!=NULL)
750 {
751 groups->ngrpnr[i]=state->natoms;
752 snew(new_egrp[i],state->natoms);
753 }
754 }
755
756 rm=0;
757 for (i=0;i<state->natoms+n;i++)
758 {
759 bRM=FALSE;
760 for(j=0;j<n;j++)
761 {
762 if(i==list[j])
763 {
764 bRM=TRUE;
765 rm++;
766 }
767 }
768
769 if(!bRM)
770 {
771 for(j=0;j<egcNR;j++)
772 {
773 if(groups->grpnr[j]!=NULL)
774 {
775 new_egrp[j][i-rm]=groups->grpnr[j][i];
776 }
777 }
778 copy_rvec(state->x[i],x_tmp[i-rm]);
779 copy_rvec(state->v[i],v_tmp[i-rm]);
780 for(j=0;j<ins_at->nr;j++)
781 {
782 if (i==ins_at->index[j])
783 ins_at->index[j]=i-rm;
784 }
785 for(j=0;j<pos_ins->pieces;j++)
786 {
787 for(k=0;k<pos_ins->nidx[j];k++)
788 {
789 if (i==pos_ins->subindex[j][k])
790 pos_ins->subindex[j][k]=i-rm;
791 }
792 }
793 }
794 }
795 sfree(state->x);
796 state->x=x_tmp;
797 sfree(state->v);
798 state->v=v_tmp;
799
800 for(i=0;i<egcNR;i++)
801 {
802 if(groups->grpnr[i]!=NULL)
803 {
804 sfree(groups->grpnr[i]);
805 groups->grpnr[i]=new_egrp[i];
806 }
807 }
808 }
809
810 int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop)
811 {
812 int i,j,m;
813 int type,natom,nmol,at,atom1=0,rm_at=0;
814 gmx_bool *bRM,bINS;
815 /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/
816 /*this routine does not live as dangerously as it seems. There is namely a check in mdrunner_membed to make
817 *sure that g_membed exits with a warning when there are molecules of the same type not in the
818 *ins_at index group. MGWolf 050710 */
819
820
821 snew(bRM,mtop->nmoltype);
822 for (i=0;i<mtop->nmoltype;i++)
823 {
824 bRM[i]=TRUE;
825 }
826
827 for (i=0;i<mtop->nmolblock;i++)
828 {
829 /*loop over molecule blocks*/
830 type =mtop->molblock[i].type;
831 natom =mtop->molblock[i].natoms_mol;
832 nmol =mtop->molblock[i].nmol;
833
834 for(j=0;j<natom*nmol && bRM[type]==TRUE;j++)
835 {
836 /*loop over atoms in the block*/
837 at=j+atom1; /*atom index = block index + offset*/
838 bINS=FALSE;
839
840 for (m=0;(m<ins_at->nr) && (bINS==FALSE);m++)
841 {
842 /*loop over atoms in insertion index group to determine if we're inserting one*/
843 if(at==ins_at->index[m])
844 {
845 bINS=TRUE;
846 }
847 }
848 bRM[type]=bINS;
849 }
850 atom1+=natom*nmol; /*update offset*/
851 if(bRM[type])
852 {
853 rm_at+=natom*nmol; /*increment bonded removal counter by # atoms in block*/
854 }
855 }
856
857 for(i=0;i<mtop->nmoltype;i++)
858 {
859 if(bRM[i])
860 {
861 for(j=0;j<F_LJ;j++)
862 {
863 mtop->moltype[i].ilist[j].nr=0;
864 }
865 for(j=F_POSRES;j<=F_VSITEN;j++)
866 {
867 mtop->moltype[i].ilist[j].nr=0;
868 }
869 }
870 }
871 sfree(bRM);
872
873 return rm_at;
874 }
875
876 void top_update(const char *topfile, char *ins, rmm_t *rm_p, gmx_mtop_t *mtop)
877 {
878 #define TEMP_FILENM "temp.top"
879 int bMolecules=0;
880 FILE *fpin,*fpout;
881 char buf[STRLEN],buf2[STRLEN],*temp;
882 int i,*nmol_rm,nmol,line;
883
884 fpin = ffopen(topfile,"r");
885 fpout = ffopen(TEMP_FILENM,"w");
886
887 snew(nmol_rm,mtop->nmoltype);
888 for(i=0;i<rm_p->nr;i++)
889 nmol_rm[rm_p->block[i]]++;
890
891 line=0;
892 while(fgets(buf,STRLEN,fpin))
893 {
894 line++;
895 if(buf[0]!=';')
896 {
897 strcpy(buf2,buf);
898 if ((temp=strchr(buf2,'\n')) != NULL)
899 temp[0]='\0';
900 ltrim(buf2);
901
902 if (buf2[0]=='[')
903 {
904 buf2[0]=' ';
905 if ((temp=strchr(buf2,'\n')) != NULL)
906 temp[0]='\0';
907 rtrim(buf2);
908 if (buf2[strlen(buf2)-1]==']')
909 {
910 buf2[strlen(buf2)-1]='\0';
911 ltrim(buf2);
912 rtrim(buf2);
913 if (gmx_strcasecmp(buf2,"molecules")==0)
914 bMolecules=1;
915 }
916 fprintf(fpout,"%s",buf);
917 } else if (bMolecules==1)
918 {
919 for(i=0;i<mtop->nmolblock;i++)
920 {
921 nmol=mtop->molblock[i].nmol;
922 sprintf(buf,"%-15s %5d\n",*(mtop->moltype[mtop->molblock[i].type].name),nmol);
923 fprintf(fpout,"%s",buf);
924 }
925 bMolecules=2;
926 } else if (bMolecules==2)
927 {
928 /* print nothing */
929 } else
930 {
931 fprintf(fpout,"%s",buf);
932 }
933 } else
934 {
935 fprintf(fpout,"%s",buf);
936 }
937 }
938
939 fclose(fpout);
940 /* use ffopen to generate backup of topinout */
941 fpout=ffopen(topfile,"w");
942 fclose(fpout);
943 rename(TEMP_FILENM,topfile);
944 #undef TEMP_FILENM
945 }
946
947 double do_md_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
948 const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
949 int nstglobalcomm,
950 gmx_vsite_t *vsite,gmx_constr_t constr,
951 int stepout,t_inputrec *ir,
952 gmx_mtop_t *top_global,
953 t_fcdata *fcd,
954 t_state *state_global,
955 t_mdatoms *mdatoms,
956 t_nrnb *nrnb,gmx_wallcycle_t wcycle,
957 gmx_edsam_t ed,t_forcerec *fr,
958 int repl_ex_nst,int repl_ex_seed,
959 real cpt_period,real max_hours,
960 const char *deviceOptions,
961 unsigned long Flags,
962 gmx_runtime_t *runtime,
963 rvec fac, rvec *r_ins, pos_ins_t *pos_ins, t_block *ins_at,
964 real xy_step, real z_step, int it_xy, int it_z)
965 {
966 gmx_mdoutf_t *outf;
967 gmx_large_int_t step,step_rel;
968 double run_time;
969 double t,t0,lam0;
970 gmx_bool bGStatEveryStep,bGStat,bNstEner,bCalcEnerPres;
971 gmx_bool bNS,bNStList,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE,
972 bFirstStep,bStateFromTPX,bInitStep,bLastStep,
973 bBornRadii,bStartingFromCpt;
974 gmx_bool bDoDHDL=FALSE;
975 gmx_bool do_ene,do_log,do_verbose,bRerunWarnNoV=TRUE,
976 bForceUpdate=FALSE,bCPT;
977 int mdof_flags;
978 gmx_bool bMasterState;
979 int force_flags,cglo_flags;
980 tensor force_vir,shake_vir,total_vir,tmp_vir,pres;
981 int i,m;
982 t_trxstatus *status;
983 rvec mu_tot;
984 t_vcm *vcm;
985 t_state *bufstate=NULL;
986 matrix *scale_tot,pcoupl_mu,M,ebox;
987 gmx_nlheur_t nlh;
988 t_trxframe rerun_fr;
989 /* gmx_repl_ex_t repl_ex=NULL;*/
990 int nchkpt=1;
991
992 gmx_localtop_t *top;
993 t_mdebin *mdebin=NULL;
994 t_state *state=NULL;
995 rvec *f_global=NULL;
996 int n_xtc=-1;
997 rvec *x_xtc=NULL;
998 gmx_enerdata_t *enerd;
999 rvec *f=NULL;
1000 gmx_global_stat_t gstat;
1001 gmx_update_t upd=NULL;
1002 t_graph *graph=NULL;
1003 globsig_t gs;
1004
1005 gmx_bool bFFscan;
1006 gmx_groups_t *groups;
1007 gmx_ekindata_t *ekind, *ekind_save;
1008 gmx_shellfc_t shellfc;
1009 int count,nconverged=0;
1010 real timestep=0;
1011 double tcount=0;
1012 gmx_bool bIonize=FALSE;
1013 gmx_bool bTCR=FALSE,bConverged=TRUE,bOK,bSumEkinhOld,bExchanged;
1014 gmx_bool bAppend;
1015 gmx_bool bResetCountersHalfMaxH=FALSE;
1016 gmx_bool bVV,bIterations,bIterate,bFirstIterate,bTemp,bPres,bTrotter;
1017 real temp0,dvdl;
1018 int a0,a1,ii;
1019 rvec *xcopy=NULL,*vcopy=NULL,*cbuf=NULL;
1020 matrix boxcopy={{0}},lastbox;
1021 real veta_save,pcurr,scalevir,tracevir;
1022 real vetanew = 0;
1023 double cycles;
1024 real last_conserved = 0;
1025 real last_ekin = 0;
1026 t_extmass MassQ;
1027 int **trotter_seq;
1028 char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
1029 int handled_stop_condition=gmx_stop_cond_none; /* compare to get_stop_condition*/
1030 gmx_iterate_t iterate;
1031 #ifdef GMX_FAHCORE
1032 /* Temporary addition for FAHCORE checkpointing */
1033 int chkpt_ret;
1034 #endif
1035
1036 /* Check for special mdrun options */
1037 bRerunMD = (Flags & MD_RERUN);
1038 bIonize = (Flags & MD_IONIZE);
1039 bFFscan = (Flags & MD_FFSCAN);
1040 bAppend = (Flags & MD_APPENDFILES);
1041 bGStatEveryStep = FALSE;
1042 if (Flags & MD_RESETCOUNTERSHALFWAY)
1043 {
1044 if (ir->nsteps > 0)
1045 {
1046 /* Signal to reset the counters half the simulation steps. */
1047 wcycle_set_reset_counters(wcycle,ir->nsteps/2);
1048 }
1049 /* Signal to reset the counters halfway the simulation time. */
1050 bResetCountersHalfMaxH = (max_hours > 0);
1051 }
1052
1053 /* md-vv uses averaged full step velocities for T-control
1054 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
1055 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
1056 bVV = EI_VV(ir->eI);
1057 if (bVV) /* to store the initial velocities while computing virial */
1058 {
1059 snew(cbuf,top_global->natoms);
1060 }
1061 /* all the iteratative cases - only if there are constraints */
1062 bIterations = ((IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
1063 bTrotter = (bVV && (IR_NPT_TROTTER(ir) || (IR_NVT_TROTTER(ir))));
1064
1065 if (bRerunMD)
1066 {
1067 /* Since we don't know if the frames read are related in any way,
1068 * rebuild the neighborlist at every step.
1069 */
1070 ir->nstlist = 1;
1071 ir->nstcalcenergy = 1;
1072 nstglobalcomm = 1;
1073 }
1074
1075 check_ir_old_tpx_versions(cr,fplog,ir,top_global);
1076
1077 nstglobalcomm = check_nstglobalcomm(fplog,cr,nstglobalcomm,ir);
1078 /*bGStatEveryStep = (nstglobalcomm == 1);*/
1079 bGStatEveryStep = FALSE;
1080
1081 if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
1082 {
1083 fprintf(fplog,
1084 "To reduce the energy communication with nstlist = -1\n"
1085 "the neighbor list validity should not be checked at every step,\n"
1086 "this means that exact integration is not guaranteed.\n"
1087 "The neighbor list validity is checked after:\n"
1088 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
1089 "In most cases this will result in exact integration.\n"
1090 "This reduces the energy communication by a factor of 2 to 3.\n"
1091 "If you want less energy communication, set nstlist > 3.\n\n");
1092 }
1093
1094 if (bRerunMD || bFFscan)
1095 {
1096 ir->nstxtcout = 0;
1097 }
1098 groups = &top_global->groups;
1099
1100 /* Initial values */
1101 init_md(fplog,cr,ir,oenv,&t,&t0,&state_global->lambda,&lam0,
1102 nrnb,top_global,&upd,
1103 nfile,fnm,&outf,&mdebin,
1104 force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
1105
1106 clear_mat(total_vir);
1107 clear_mat(pres);
1108 /* Energy terms and groups */
1109 snew(enerd,1);
1110 init_enerdata(top_global->groups.grps[egcENER].nr,ir->n_flambda,enerd);
1111 if (DOMAINDECOMP(cr))
1112 {
1113 f = NULL;
1114 }
1115 else
1116 {
1117 snew(f,top_global->natoms);
1118 }
1119
1120 /* Kinetic energy data */
1121 snew(ekind,1);
1122 init_ekindata(fplog,top_global,&(ir->opts),ekind);
1123 /* needed for iteration of constraints */
1124 snew(ekind_save,1);
1125 init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
1126 /* Copy the cos acceleration to the groups struct */
1127 ekind->cosacc.cos_accel = ir->cos_accel;
1128
1129 gstat = global_stat_init(ir);
1130 debug_gmx();
1131
1132 /* Check for polarizable models and flexible constraints */
1133 shellfc = init_shell_flexcon(fplog,
1134 top_global,n_flexible_constraints(constr),
1135 (ir->bContinuation ||
1136 (DOMAINDECOMP(cr) && !MASTER(cr))) ?
1137 NULL : state_global->x);
1138
1139 /* if (DEFORM(*ir))
1140 {
1141 #ifdef GMX_THREADS
1142 tMPI_Thread_mutex_lock(&deform_init_box_mutex);
1143 #endif
1144 set_deform_reference_box(upd,
1145 deform_init_init_step_tpx,
1146 deform_init_box_tpx);
1147 #ifdef GMX_THREADS
1148 tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
1149 #endif
1150 }*/
1151
1152 /* {
1153 double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
1154 if ((io > 2000) && MASTER(cr))
1155 fprintf(stderr,
1156 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
1157 io);
1158 }*/
1159
1160 if (DOMAINDECOMP(cr)) {
1161 top = dd_init_local_top(top_global);
1162
1163 snew(state,1);
1164 dd_init_local_state(cr->dd,state_global,state);
1165
1166 if (DDMASTER(cr->dd) && ir->nstfout) {
1167 snew(f_global,state_global->natoms);
1168 }
1169 } else {
1170 if (PAR(cr)) {
1171 /* Initialize the particle decomposition and split the topology */
1172 top = split_system(fplog,top_global,ir,cr);
1173
1174 pd_cg_range(cr,&fr->cg0,&fr->hcg);
1175 pd_at_range(cr,&a0,&a1);
1176 } else {
1177 top = gmx_mtop_generate_local_top(top_global,ir);
1178
1179 a0 = 0;
1180 a1 = top_global->natoms;
1181 }
1182
1183 state = partdec_init_local_state(cr,state_global);
1184 f_global = f;
1185
1186 atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
1187
1188 if (vsite) {
1189 set_vsite_top(vsite,top,mdatoms,cr);
1190 }
1191
1192 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols) {
1193 graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
1194 }
1195
1196 if (shellfc) {
1197 make_local_shells(cr,mdatoms,shellfc);
1198 }
1199
1200 if (ir->pull && PAR(cr)) {
1201 dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
1202 }
1203 }
1204
1205 if (DOMAINDECOMP(cr))
1206 {
1207 /* Distribute the charge groups over the nodes from the master node */
1208 dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
1209 state_global,top_global,ir,
1210 state,&f,mdatoms,top,fr,
1211 vsite,shellfc,constr,
1212 nrnb,wcycle,FALSE);
1213 }
1214
1215 update_mdatoms(mdatoms,state->lambda);
1216
1217 if (MASTER(cr))
1218 {
1219 if (opt2bSet("-cpi",nfile,fnm))
1220 {
1221 /* Update mdebin with energy history if appending to output files */
1222 if ( Flags & MD_APPENDFILES )
1223 {
1224 restore_energyhistory_from_state(mdebin,&state_global->enerhist);
1225 }
1226 else
1227 {
1228 /* We might have read an energy history from checkpoint,
1229 * free the allocated memory and reset the counts.
1230 */
1231 done_energyhistory(&state_global->enerhist);
1232 init_energyhistory(&state_global->enerhist);
1233 }
1234 }
1235 /* Set the initial energy history in state by updating once */
1236 update_energyhistory(&state_global->enerhist,mdebin);
1237 }
1238
1239 if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG)) {
1240 /* Set the random state if we read a checkpoint file */
1241 set_stochd_state(upd,state);
1242 }
1243
1244 /* Initialize constraints */
1245 if (constr) {
1246 if (!DOMAINDECOMP(cr))
1247 set_constraints(constr,top,ir,mdatoms,cr);
1248 }
1249
1250 /* Check whether we have to GCT stuff */
1251 /* bTCR = ftp2bSet(efGCT,nfile,fnm);
1252 if (bTCR) {
1253 if (MASTER(cr)) {
1254 fprintf(stderr,"Will do General Coupling Theory!\n");
1255 }
1256 gnx = top_global->mols.nr;
1257 snew(grpindex,gnx);
1258 for(i=0; (i<gnx); i++) {
1259 grpindex[i] = i;
1260 }
1261 }*/
1262
1263 /* if (repl_ex_nst > 0 && MASTER(cr))
1264 repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
1265 repl_ex_nst,repl_ex_seed);*/
1266
1267 if (!ir->bContinuation && !bRerunMD)
1268 {
1269 if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
1270 {
1271 /* Set the velocities of frozen particles to zero */
1272 for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
1273 {
1274 for(m=0; m<DIM; m++)
1275 {
1276 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
1277 {
1278 state->v[i][m] = 0;
1279 }
1280 }
1281 }
1282 }
1283
1284 if (constr)
1285 {
1286 /* Constrain the initial coordinates and velocities */
1287 do_constrain_first(fplog,constr,ir,mdatoms,state,f,
1288 graph,cr,nrnb,fr,top,shake_vir);
1289 }
1290 if (vsite)
1291 {
1292 /* Construct the virtual sites for the initial configuration */
1293 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
1294 top->idef.iparams,top->idef.il,
1295 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1296 }
1297 }
1298
1299 debug_gmx();
1300
1301 /* I'm assuming we need global communication the first time! MRS */
1302 cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
1303 | (bVV ? CGLO_PRESSURE:0)
1304 | (bVV ? CGLO_CONSTRAINT:0)
1305 | (bRerunMD ? CGLO_RERUNMD:0)
1306 | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
1307
1308 bSumEkinhOld = FALSE;
1309 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1310 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1311 constr,NULL,FALSE,state->box,
1312 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
1313 if (ir->eI == eiVVAK) {
1314 /* a second call to get the half step temperature initialized as well */
1315 /* we do the same call as above, but turn the pressure off -- internally, this
1316 is recognized as a velocity verlet half-step kinetic energy calculation.
1317 This minimized excess variables, but perhaps loses some logic?*/
1318
1319 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1320 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1321 constr,NULL,FALSE,state->box,
1322 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1323 cglo_flags &~ CGLO_PRESSURE);
1324 }
1325
1326 /* Calculate the initial half step temperature, and save the ekinh_old */
1327 if (!(Flags & MD_STARTFROMCPT))
1328 {
1329 for(i=0; (i<ir->opts.ngtc); i++)
1330 {
1331 copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
1332 }
1333 }
1334 if (ir->eI != eiVV)
1335 {
1336 enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
1337 and there is no previous step */
1338 }
1339 temp0 = enerd->term[F_TEMP];
1340
1341 /* if using an iterative algorithm, we need to create a working directory for the state. */
1342 if (bIterations)
1343 {
1344 bufstate = init_bufstate(state);
1345 }
1346 if (bFFscan)
1347 {
1348 snew(xcopy,state->natoms);
1349 snew(vcopy,state->natoms);
1350 copy_rvecn(state->x,xcopy,0,state->natoms);
1351 copy_rvecn(state->v,vcopy,0,state->natoms);
1352 copy_mat(state->box,boxcopy);
1353 }
1354
1355 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
1356 temperature control */
1357 trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
1358
1359 if (MASTER(cr))
1360 {
1361 if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
1362 {
1363 fprintf(fplog,
1364 "RMS relative constraint deviation after constraining: %.2e\n",
1365 constr_rmsd(constr,FALSE));
1366 }
1367 fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
1368 if (bRerunMD)
1369 {
1370 fprintf(stderr,"starting md rerun '%s', reading coordinates from"
1371 " input trajectory '%s'\n\n",
1372 *(top_global->name),opt2fn("-rerun",nfile,fnm));
1373 if (bVerbose)
1374 {
1375 fprintf(stderr,"Calculated time to finish depends on nsteps from "
1376 "run input file,\nwhich may not correspond to the time "
1377 "needed to process input trajectory.\n\n");
1378 }
1379 }
1380 else
1381 {
1382 char tbuf[20];
1383 fprintf(stderr,"starting mdrun '%s'\n",
1384 *(top_global->name));
1385 if (ir->nsteps >= 0)
1386 {
1387 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
1388 }
1389 else
1390 {
1391 sprintf(tbuf,"%s","infinite");
1392 }
1393 if (ir->init_step > 0)
1394 {
1395 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
1396 gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
1397 gmx_step_str(ir->init_step,sbuf2),
1398 ir->init_step*ir->delta_t);
1399 }
1400 else
1401 {
1402 fprintf(stderr,"%s steps, %s ps.\n",
1403 gmx_step_str(ir->nsteps,sbuf),tbuf);
1404 }
1405 }
1406 fprintf(fplog,"\n");
1407 }
1408
1409 /* Set and write start time */
1410 runtime_start(runtime);
1411 print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
1412 wallcycle_start(wcycle,ewcRUN);
1413 if (fplog)
1414 fprintf(fplog,"\n");
1415
1416 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
1417 /*#ifdef GMX_FAHCORE
1418 chkpt_ret=fcCheckPointParallel( cr->nodeid,
1419 NULL,0);
1420 if ( chkpt_ret == 0 )
1421 gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
1422 #endif*/
1423
1424 debug_gmx();
1425 /***********************************************************
1426 *
1427 * Loop over MD steps
1428 *
1429 ************************************************************/
1430
1431 /* if rerunMD then read coordinates and velocities from input trajectory */
1432 if (bRerunMD)
1433 {
1434 if (getenv("GMX_FORCE_UPDATE"))
1435 {
1436 bForceUpdate = TRUE;
1437 }
1438
1439 bNotLastFrame = read_first_frame(oenv,&status,
1440 opt2fn("-rerun",nfile,fnm),
1441 &rerun_fr,TRX_NEED_X | TRX_READ_V);
1442 if (rerun_fr.natoms != top_global->natoms)
1443 {
1444 gmx_fatal(FARGS,
1445 "Number of atoms in trajectory (%d) does not match the "
1446 "run input file (%d)\n",
1447 rerun_fr.natoms,top_global->natoms);
1448 }
1449 if (ir->ePBC != epbcNONE)
1450 {
1451 if (!rerun_fr.bBox)
1452 {
1453 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);
1454 }
1455 if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
1456 {
1457 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
1458 }
1459
1460 /* Set the shift vectors.
1461 * Necessary here when have a static box different from the tpr box.
1462 */
1463 calc_shifts(rerun_fr.box,fr->shift_vec);
1464 }
1465 }
1466
1467 /* loop over MD steps or if rerunMD to end of input trajectory */
1468 bFirstStep = TRUE;
1469 /* Skip the first Nose-Hoover integration when we get the state from tpx */
1470 bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
1471 bInitStep = bFirstStep && (bStateFromTPX || bVV);
1472 bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
1473 bLastStep = FALSE;
1474 bSumEkinhOld = FALSE;
1475 bExchanged = FALSE;
1476
1477 init_global_signals(&gs,cr,ir,repl_ex_nst);
1478
1479 step = ir->init_step;
1480 step_rel = 0;
1481
1482 if (ir->nstlist == -1)
1483 {
1484 init_nlistheuristics(&nlh,bGStatEveryStep,step);
1485 }
1486
1487 bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps));
1488 while (!bLastStep || (bRerunMD && bNotLastFrame)) {
1489
1490 wallcycle_start(wcycle,ewcSTEP);
1491
1492 GMX_MPE_LOG(ev_timestep1);
1493
1494 if (bRerunMD) {
1495 if (rerun_fr.bStep) {
1496 step = rerun_fr.step;
1497 step_rel = step - ir->init_step;
1498 }
1499 if (rerun_fr.bTime) {
1500 t = rerun_fr.time;
1501 }
1502 else
1503 {
1504 t = step;
1505 }
1506 }
1507 else
1508 {
1509 bLastStep = (step_rel == ir->nsteps);
1510 t = t0 + step*ir->delta_t;
1511 }
1512
1513 if (ir->efep != efepNO)
1514 {
1515 if (bRerunMD && rerun_fr.bLambda && (ir->delta_lambda!=0))
1516 {
1517 state_global->lambda = rerun_fr.lambda;
1518 }
1519 else
1520 {
1521 state_global->lambda = lam0 + step*ir->delta_lambda;
1522 }
1523 state->lambda = state_global->lambda;
1524 bDoDHDL = do_per_step(step,ir->nstdhdl);
1525 }
1526
1527 if (bSimAnn)
1528 {
1529 update_annealing_target_temp(&(ir->opts),t);
1530 }
1531
1532 if (bRerunMD)
1533 {
1534 if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
1535 {
1536 for(i=0; i<state_global->natoms; i++)
1537 {
1538 copy_rvec(rerun_fr.x[i],state_global->x[i]);
1539 }
1540 if (rerun_fr.bV)
1541 {
1542 for(i=0; i<state_global->natoms; i++)
1543 {
1544 copy_rvec(rerun_fr.v[i],state_global->v[i]);
1545 }
1546 }
1547 else
1548 {
1549 for(i=0; i<state_global->natoms; i++)
1550 {
1551 clear_rvec(state_global->v[i]);
1552 }
1553 if (bRerunWarnNoV)
1554 {
1555 fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
1556 " Ekin, temperature and pressure are incorrect,\n"
1557 " the virial will be incorrect when constraints are present.\n"
1558 "\n");
1559 bRerunWarnNoV = FALSE;
1560 }
1561 }
1562 }
1563 copy_mat(rerun_fr.box,state_global->box);
1564 copy_mat(state_global->box,state->box);
1565
1566 if (vsite && (Flags & MD_RERUN_VSITE))
1567 {
1568 if (DOMAINDECOMP(cr))
1569 {
1570 gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
1571 }
1572 if (graph)
1573 {
1574 /* Following is necessary because the graph may get out of sync
1575 * with the coordinates if we only have every N'th coordinate set
1576 */
1577 mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
1578 shift_self(graph,state->box,state->x);
1579 }
1580 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
1581 top->idef.iparams,top->idef.il,
1582 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1583 if (graph)
1584 {
1585 unshift_self(graph,state->box,state->x);
1586 }
1587 }
1588 }
1589
1590 /* Stop Center of Mass motion */
1591 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
1592
1593 /* Copy back starting coordinates in case we're doing a forcefield scan */
1594 if (bFFscan)
1595 {
1596 for(ii=0; (ii<state->natoms); ii++)
1597 {
1598 copy_rvec(xcopy[ii],state->x[ii]);
1599 copy_rvec(vcopy[ii],state->v[ii]);
1600 }
1601 copy_mat(boxcopy,state->box);
1602 }
1603
1604 if (bRerunMD)
1605 {
1606 /* for rerun MD always do Neighbour Searching */
1607 bNS = (bFirstStep || ir->nstlist != 0);
1608 bNStList = bNS;
1609 }
1610 else
1611 {
1612 /* Determine whether or not to do Neighbour Searching and LR */
1613 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
1614
1615 bNS = (bFirstStep || bExchanged || bNStList ||
1616 (ir->nstlist == -1 && nlh.nabnsb > 0));
1617
1618 if (bNS && ir->nstlist == -1)
1619 {
1620 set_nlistheuristics(&nlh,bFirstStep || bExchanged,step);
1621 }
1622 }
1623
1624 /* < 0 means stop at next step, > 0 means stop at next NS step */
1625 if ( (gs.set[eglsSTOPCOND] < 0 ) ||
1626 ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
1627 {
1628 bLastStep = TRUE;
1629 }
1630
1631 /* Determine whether or not to update the Born radii if doing GB */
1632 bBornRadii=bFirstStep;
1633 if (ir->implicit_solvent && (step % ir->nstgbradii==0))
1634 {
1635 bBornRadii=TRUE;
1636 }
1637
1638 do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
1639 do_verbose = bVerbose &&
1640 (step % stepout == 0 || bFirstStep || bLastStep);
1641
1642 if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
1643 {
1644 if (bRerunMD)
1645 {
1646 bMasterState = TRUE;
1647 }
1648 else
1649 {
1650 bMasterState = FALSE;
1651 /* Correct the new box if it is too skewed */
1652 if (DYNAMIC_BOX(*ir))
1653 {
1654 if (correct_box(fplog,step,state->box,graph))
1655 {
1656 bMasterState = TRUE;
1657 }
1658 }
1659 if (DOMAINDECOMP(cr) && bMasterState)
1660 {
1661 dd_collect_state(cr->dd,state,state_global);
1662 }
1663 }
1664
1665 if (DOMAINDECOMP(cr))
1666 {
1667 /* Repartition the domain decomposition */
1668 wallcycle_start(wcycle,ewcDOMDEC);
1669 dd_partition_system(fplog,step,cr,
1670 bMasterState,nstglobalcomm,
1671 state_global,top_global,ir,
1672 state,&f,mdatoms,top,fr,
1673 vsite,shellfc,constr,
1674 nrnb,wcycle,do_verbose);
1675 wallcycle_stop(wcycle,ewcDOMDEC);
1676 /* If using an iterative integrator, reallocate space to match the decomposition */
1677 }
1678 }
1679
1680 if (MASTER(cr) && do_log && !bFFscan)
1681 {
1682 print_ebin_header(fplog,step,t,state->lambda);
1683 }
1684
1685 if (ir->efep != efepNO)
1686 {
1687 update_mdatoms(mdatoms,state->lambda);
1688 }
1689
1690 if (bRerunMD && rerun_fr.bV)
1691 {
1692
1693 /* We need the kinetic energy at minus the half step for determining
1694 * the full step kinetic energy and possibly for T-coupling.*/
1695 /* This may not be quite working correctly yet . . . . */
1696 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1697 wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
1698 constr,NULL,FALSE,state->box,
1699 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1700 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
1701 }
1702 clear_mat(force_vir);
1703
1704 /* Ionize the atoms if necessary */
1705 /* if (bIonize)
1706 {
1707 ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
1708 mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
1709 }*/
1710
1711 /* Update force field in ffscan program */
1712 /* if (bFFscan)
1713 {
1714 if (update_forcefield(fplog,
1715 nfile,fnm,fr,
1716 mdatoms->nr,state->x,state->box)) {
1717 if (gmx_parallel_env_initialized())
1718 {
1719 gmx_finalize();
1720 }
1721 exit(0);
1722 }
1723 }*/
1724
1725 GMX_MPE_LOG(ev_timestep2);
1726
1727 /* We write a checkpoint at this MD step when:
1728 * either at an NS step when we signalled through gs,
1729 * or at the last step (but not when we do not want confout),
1730 * but never at the first step or with rerun.
1731 */
1732 /* bCPT = (((gs.set[eglsCHKPT] && bNS) ||
1733 (bLastStep && (Flags & MD_CONFOUT))) &&
1734 step > ir->init_step && !bRerunMD);
1735 if (bCPT)
1736 {
1737 gs.set[eglsCHKPT] = 0;
1738 }*/
1739
1740 /* Determine the energy and pressure:
1741 * at nstcalcenergy steps and at energy output steps (set below).
1742 */
1743 bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
1744 bCalcEnerPres = bNstEner;
1745
1746 /* Do we need global communication ? */
1747 bGStat = (bCalcEnerPres || bStopCM ||
1748 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
1749
1750 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
1751
1752 if (do_ene || do_log)
1753 {
1754 bCalcEnerPres = TRUE;
1755 bGStat = TRUE;
1756 }
1757
1758 /* these CGLO_ options remain the same throughout the iteration */
1759 cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
1760 (bStopCM ? CGLO_STOPCM : 0) |
1761 (bGStat ? CGLO_GSTAT : 0)
1762 );
1763
1764 force_flags = (GMX_FORCE_STATECHANGED |
1765 ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
1766 GMX_FORCE_ALLFORCES |
1767 (bNStList ? GMX_FORCE_DOLR : 0) |
1768 GMX_FORCE_SEPLRF |
1769 (bCalcEnerPres ? GMX_FORCE_VIRIAL : 0) |
1770 (bDoDHDL ? GMX_FORCE_DHDL : 0)
1771 );
1772
1773 if (shellfc)
1774 {
1775 /* Now is the time to relax the shells */
1776 count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
1777 ir,bNS,force_flags,
1778 bStopCM,top,top_global,
1779 constr,enerd,fcd,
1780 state,f,force_vir,mdatoms,
1781 nrnb,wcycle,graph,groups,
1782 shellfc,fr,bBornRadii,t,mu_tot,
1783 state->natoms,&bConverged,vsite,
1784 outf->fp_field);
1785 tcount+=count;
1786
1787 if (bConverged)
1788 {
1789 nconverged++;
1790 }
1791 }
1792 else
1793 {
1794 /* The coordinates (x) are shifted (to get whole molecules)
1795 * in do_force.
1796 * This is parallellized as well, and does communication too.
1797 * Check comments in sim_util.c
1798 */
1799
1800 do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
1801 state->box,state->x,&state->hist,
1802 f,force_vir,mdatoms,enerd,fcd,
1803 state->lambda,graph,
1804 fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
1805 (bNS ? GMX_FORCE_NS : 0) | force_flags);
1806 }
1807
1808 GMX_BARRIER(cr->mpi_comm_mygroup);
1809
1810 /* if (bTCR)
1811 {
1812 mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
1813 mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
1814 }
1815
1816 if (bTCR && bFirstStep)
1817 {
1818 tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
1819 fprintf(fplog,"Done init_coupling\n");
1820 fflush(fplog);
1821 }*/
1822
1823 /* ############### START FIRST UPDATE HALF-STEP ############### */
1824
1825 if (bVV && !bStartingFromCpt && !bRerunMD)
1826 {
1827 if (ir->eI == eiVV)
1828 {
1829 if (bInitStep)
1830 {
1831 /* if using velocity verlet with full time step Ekin,
1832 * take the first half step only to compute the
1833 * virial for the first step. From there,
1834 * revert back to the initial coordinates
1835 * so that the input is actually the initial step.
1836 */
1837 copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
1838 }
1839
1840 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1841 if (!bInitStep)
1842 {
1843 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
1844 }
1845
1846 if (ir->eI == eiVVAK)
1847 {
1848 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
1849 }
1850
1851 update_coords(fplog,step,ir,mdatoms,state,
1852 f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1853 ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
1854 cr,nrnb,constr,&top->idef);
1855
1856 if (bIterations)
1857 {
1858 gmx_iterate_init(&iterate,bIterations && !bInitStep);
1859 }
1860 /* for iterations, we save these vectors, as we will be self-consistently iterating
1861 the calculations */
1862 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1863
1864 /* save the state */
1865 if (bIterations && iterate.bIterate) {
1866 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1867 }
1868 }
1869
1870 bFirstIterate = TRUE;
1871 while (bFirstIterate || (bIterations && iterate.bIterate))
1872 {
1873 if (bIterations && iterate.bIterate)
1874 {
1875 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1876 if (bFirstIterate && bTrotter)
1877 {
1878 /* The first time through, we need a decent first estimate
1879 of veta(t+dt) to compute the constraints. Do
1880 this by computing the box volume part of the
1881 trotter integration at this time. Nothing else
1882 should be changed by this routine here. If
1883 !(first time), we start with the previous value
1884 of veta. */
1885
1886 veta_save = state->veta;
1887 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
1888 vetanew = state->veta;
1889 state->veta = veta_save;
1890 }
1891 }
1892
1893 bOK = TRUE;
1894 if ( !bRerunMD || rerun_fr.bV || bForceUpdate) { /* Why is rerun_fr.bV here? Unclear. */
1895 dvdl = 0;
1896
1897 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1898 &top->idef,shake_vir,NULL,
1899 cr,nrnb,wcycle,upd,constr,
1900 bInitStep,TRUE,bCalcEnerPres,vetanew);
1901
1902 if (!bOK && !bFFscan)
1903 {
1904 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1905 }
1906
1907 }
1908 else if (graph)
1909 { /* Need to unshift here if a do_force has been
1910 called in the previous step */
1911 unshift_self(graph,state->box,state->x);
1912 }
1913
1914
1915 if (bVV) {
1916 /* if VV, compute the pressure and constraints */
1917 /* if VV2, the pressure and constraints only if using pressure control.*/
1918 bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir));
1919 bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));
1920 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1921 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1922 constr,NULL,FALSE,state->box,
1923 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1924 cglo_flags
1925 | CGLO_ENERGY
1926 | (bTemp ? CGLO_TEMPERATURE:0)
1927 | (bPres ? CGLO_PRESSURE : 0)
1928 | (bPres ? CGLO_CONSTRAINT : 0)
1929 | (iterate.bIterate ? CGLO_ITERATE : 0)
1930 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1931 | CGLO_SCALEEKIN
1932 );
1933 }
1934 /* explanation of above:
1935 a) We compute Ekin at the full time step
1936 if 1) we are using the AveVel Ekin, and it's not the
1937 initial step, or 2) if we are using AveEkin, but need the full
1938 time step kinetic energy for the pressure.
1939 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
1940 EkinAveVel because it's needed for the pressure */
1941
1942 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1943 if (bVV && !bInitStep)
1944 {
1945 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ2);
1946 }
1947
1948 if (bIterations &&
1949 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
1950 state->veta,&vetanew))
1951 {
1952 break;
1953 }
1954 bFirstIterate = FALSE;
1955 }
1956
1957 if (bTrotter && !bInitStep) {
1958 copy_mat(shake_vir,state->svir_prev);
1959 copy_mat(force_vir,state->fvir_prev);
1960 if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
1961 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1962 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
1963 enerd->term[F_EKIN] = trace(ekind->ekin);
1964 }
1965 }
1966 /* if it's the initial step, we performed this first step just to get the constraint virial */
1967 if (bInitStep && ir->eI==eiVV) {
1968 copy_rvecn(cbuf,state->v,0,state->natoms);
1969 }
1970
1971 if (fr->bSepDVDL && fplog && do_log)
1972 {
1973 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
1974 }
1975 enerd->term[F_DHDL_CON] += dvdl;
1976
1977 GMX_MPE_LOG(ev_timestep1);
1978
1979 }
1980
1981 /* MRS -- now done iterating -- compute the conserved quantity */
1982 if (bVV) {
1983 last_conserved = 0;
1984 if (IR_NVT_TROTTER(ir) || IR_NPT_TROTTER(ir))
1985 {
1986 last_conserved =
1987 NPT_energy(ir,state,&MassQ);
1988 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1989 {
1990 last_conserved -= enerd->term[F_DISPCORR];
1991 }
1992 }
1993 if (ir->eI==eiVV) {
1994 last_ekin = enerd->term[F_EKIN]; /* does this get preserved through checkpointing? */
1995 }
1996 }
1997
1998 /* ######## END FIRST UPDATE STEP ############## */
1999 /* ######## If doing VV, we now have v(dt) ###### */
2000
2001 /* ################## START TRAJECTORY OUTPUT ################# */
2002
2003 /* Now we have the energies and forces corresponding to the
2004 * coordinates at time t. We must output all of this before
2005 * the update.
2006 * for RerunMD t is read from input trajectory
2007 */
2008 GMX_MPE_LOG(ev_output_start);
2009
2010 mdof_flags = 0;
2011 if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
2012 if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
2013 if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
2014 if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
2015 /* if (bCPT) { mdof_flags |= MDOF_CPT; };*/
2016
2017 #ifdef GMX_FAHCORE
2018 if (MASTER(cr))
2019 fcReportProgress( ir->nsteps, step );
2020
2021 if (bLastStep)
2022 {
2023 /* Enforce writing positions and velocities at end of run */
2024 mdof_flags |= (MDOF_X | MDOF_V);
2025 }
2026 /* sync bCPT and fc record-keeping */
2027 /* if (bCPT && MASTER(cr))
2028 fcRequestCheckPoint();*/
2029 #endif
2030
2031 if (mdof_flags != 0)
2032 {
2033 wallcycle_start(wcycle,ewcTRAJ);
2034 /* if (bCPT)
2035 {
2036 if (state->flags & (1<<estLD_RNG))
2037 {
2038 get_stochd_state(upd,state);
2039 }
2040 if (MASTER(cr))
2041 {
2042 if (bSumEkinhOld)
2043 {
2044 state_global->ekinstate.bUpToDate = FALSE;
2045 }
2046 else
2047 {
2048 update_ekinstate(&state_global->ekinstate,ekind);
2049 state_global->ekinstate.bUpToDate = TRUE;
2050 }
2051 update_energyhistory(&state_global->enerhist,mdebin);
2052 }
2053 }*/
2054 write_traj(fplog,cr,outf,mdof_flags,top_global,
2055 step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
2056 /* if (bCPT)
2057 {
2058 nchkpt++;
2059 bCPT = FALSE;
2060 }*/
2061 debug_gmx();
2062 if (bLastStep && step_rel == ir->nsteps &&
2063 (Flags & MD_CONFOUT) && MASTER(cr) &&
2064 !bRerunMD && !bFFscan)
2065 {
2066 /* x and v have been collected in write_traj,
2067 * because a checkpoint file will always be written
2068 * at the last step.
2069 */
2070 fprintf(stderr,"\nWriting final coordinates.\n");
2071 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
2072 DOMAINDECOMP(cr))
2073 {
2074 /* Make molecules whole only for confout writing */
2075 do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
2076 }
2077 /* write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
2078 *top_global->name,top_global,
2079 state_global->x,state_global->v,
2080 ir->ePBC,state->box);*/
2081 debug_gmx();
2082 }
2083 wallcycle_stop(wcycle,ewcTRAJ);
2084 }
2085 GMX_MPE_LOG(ev_output_finish);
2086
2087 /* kludge -- virial is lost with restart for NPT control. Must restart */
2088 if (bStartingFromCpt && bVV)
2089 {
2090 copy_mat(state->svir_prev,shake_vir);
2091 copy_mat(state->fvir_prev,force_vir);
2092 }
2093 /* ################## END TRAJECTORY OUTPUT ################ */
2094
2095 /* Determine the pressure:
2096 * always when we want exact averages in the energy file,
2097 * at ns steps when we have pressure coupling,
2098 * otherwise only at energy output steps (set below).
2099 */
2100
2101 bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
2102 bCalcEnerPres = bNstEner;
2103
2104 /* Do we need global communication ? */
2105 bGStat = (bGStatEveryStep || bStopCM || bNS ||
2106 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
2107
2108 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
2109
2110 if (do_ene || do_log)
2111 {
2112 bCalcEnerPres = TRUE;
2113 bGStat = TRUE;
2114 }
2115
2116 /* Determine the wallclock run time up till now */
2117 run_time = gmx_gettime() - (double)runtime->real;
2118
2119 /* Check whether everything is still allright */
2120 if (((int)gmx_get_stop_condition() > handled_stop_condition)
2121 #ifdef GMX_THREADS
2122 && MASTER(cr)
2123 #endif
2124 )
2125 {
2126 /* this is just make gs.sig compatible with the hack
2127 of sending signals around by MPI_Reduce with together with
2128 other floats */
2129 if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
2130 gs.sig[eglsSTOPCOND]=1;
2131 if ( gmx_get_stop_condition() == gmx_stop_cond_next )
2132 gs.sig[eglsSTOPCOND]=-1;
2133 /* < 0 means stop at next step, > 0 means stop at next NS step */
2134 if (fplog)
2135 {
2136 fprintf(fplog,
2137 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2138 gmx_get_signal_name(),
2139 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
2140 fflush(fplog);
2141 }
2142 fprintf(stderr,
2143 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2144 gmx_get_signal_name(),
2145 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
2146 fflush(stderr);
2147 handled_stop_condition=(int)gmx_get_stop_condition();
2148 }
2149 else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
2150 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
2151 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
2152 {
2153 /* Signal to terminate the run */
2154 gs.sig[eglsSTOPCOND] = 1;
2155 if (fplog)
2156 {
2157 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
2158 }
2159 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
2160 }
2161
2162 if (bResetCountersHalfMaxH && MASTER(cr) &&
2163 run_time > max_hours*60.0*60.0*0.495)
2164 {
2165 gs.sig[eglsRESETCOUNTERS] = 1;
2166 }
2167
2168 if (ir->nstlist == -1 && !bRerunMD)
2169 {
2170 /* When bGStatEveryStep=FALSE, global_stat is only called
2171 * when we check the atom displacements, not at NS steps.
2172 * This means that also the bonded interaction count check is not
2173 * performed immediately after NS. Therefore a few MD steps could
2174 * be performed with missing interactions.
2175 * But wrong energies are never written to file,
2176 * since energies are only written after global_stat
2177 * has been called.
2178 */
2179 if (step >= nlh.step_nscheck)
2180 {
2181 nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
2182 nlh.scale_tot,state->x);
2183 }
2184 else
2185 {
2186 /* This is not necessarily true,
2187 * but step_nscheck is determined quite conservatively.
2188 */
2189 nlh.nabnsb = 0;
2190 }
2191 }
2192
2193 /* In parallel we only have to check for checkpointing in steps
2194 * where we do global communication,
2195 * otherwise the other nodes don't know.
2196 */
2197 if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
2198 cpt_period >= 0 &&
2199 (cpt_period == 0 ||
2200 run_time >= nchkpt*cpt_period*60.0)) &&
2201 gs.set[eglsCHKPT] == 0)
2202 {
2203 gs.sig[eglsCHKPT] = 1;
2204 }
2205
2206 if (bIterations)
2207 {
2208 gmx_iterate_init(&iterate,bIterations);
2209 }
2210
2211 /* for iterations, we save these vectors, as we will be redoing the calculations */
2212 if (bIterations && iterate.bIterate)
2213 {
2214 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
2215 }
2216 bFirstIterate = TRUE;
2217 while (bFirstIterate || (bIterations && iterate.bIterate))
2218 {
2219 /* We now restore these vectors to redo the calculation with improved extended variables */
2220 if (bIterations)
2221 {
2222 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
2223 }
2224
2225 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
2226 so scroll down for that logic */
2227
2228 /* ######### START SECOND UPDATE STEP ################# */
2229 GMX_MPE_LOG(ev_update_start);
2230 bOK = TRUE;
2231 if (!bRerunMD || rerun_fr.bV || bForceUpdate)
2232 {
2233 wallcycle_start(wcycle,ewcUPDATE);
2234 dvdl = 0;
2235 /* Box is changed in update() when we do pressure coupling,
2236 * but we should still use the old box for energy corrections and when
2237 * writing it to the energy file, so it matches the trajectory files for
2238 * the same timestep above. Make a copy in a separate array.
2239 */
2240 copy_mat(state->box,lastbox);
2241 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
2242 if (bTrotter)
2243 {
2244 if (bIterations && iterate.bIterate)
2245 {
2246 if (bFirstIterate)
2247 {
2248 scalevir = 1;
2249 }
2250 else
2251 {
2252 /* we use a new value of scalevir to converge the iterations faster */
2253 scalevir = tracevir/trace(shake_vir);
2254 }
2255 msmul(shake_vir,scalevir,shake_vir);
2256 m_add(force_vir,shake_vir,total_vir);
2257 clear_mat(shake_vir);
2258 }
2259 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ3);
2260 }
2261 /* We can only do Berendsen coupling after we have summed
2262 * the kinetic energy or virial. Since the happens
2263 * in global_state after update, we should only do it at
2264 * step % nstlist = 1 with bGStatEveryStep=FALSE.
2265 */
2266
2267 if (ir->eI != eiVVAK)
2268 {
2269 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
2270 }
2271 update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
2272 upd,bInitStep);
2273
2274 if (bVV)
2275 {
2276 /* velocity half-step update */
2277 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2278 ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,cr,nrnb,constr,&top->idef);
2279 }
2280
2281 /* Above, initialize just copies ekinh into ekin,
2282 * it doesn't copy position (for VV),
2283 * and entire integrator for MD.
2284 */
2285
2286 if (ir->eI==eiVVAK)
2287 {
2288 copy_rvecn(state->x,cbuf,0,state->natoms);
2289 }
2290
2291 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2292 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
2293 wallcycle_stop(wcycle,ewcUPDATE);
2294
2295 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2296 &top->idef,shake_vir,force_vir,
2297 cr,nrnb,wcycle,upd,constr,
2298 bInitStep,FALSE,bCalcEnerPres,state->veta);
2299
2300 if (ir->eI==eiVVAK)
2301 {
2302 /* erase F_EKIN and F_TEMP here? */
2303 /* just compute the kinetic energy at the half step to perform a trotter step */
2304 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2305 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2306 constr,NULL,FALSE,lastbox,
2307 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2308 cglo_flags | CGLO_TEMPERATURE | CGLO_CONSTRAINT
2309 );
2310 wallcycle_start(wcycle,ewcUPDATE);
2311 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ, trotter_seq,ettTSEQ4);
2312 /* now we know the scaling, we can compute the positions again again */
2313 copy_rvecn(cbuf,state->x,0,state->natoms);
2314
2315 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2316 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
2317 wallcycle_stop(wcycle,ewcUPDATE);
2318
2319 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
2320 /* are the small terms in the shake_vir here due
2321 * to numerical errors, or are they important
2322 * physically? I'm thinking they are just errors, but not completely sure.
2323 * For now, will call without actually constraining, constr=NULL*/
2324 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2325 &top->idef,tmp_vir,force_vir,
2326 cr,nrnb,wcycle,upd,NULL,
2327 bInitStep,FALSE,bCalcEnerPres,state->veta);
2328 }
2329 if (!bOK && !bFFscan)
2330 {
2331 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2332 }
2333
2334 if (fr->bSepDVDL && fplog && do_log)
2335 {
2336 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
2337 }
2338 enerd->term[F_DHDL_CON] += dvdl;
2339 }
2340 else if (graph)
2341 {
2342 /* Need to unshift here */
2343 unshift_self(graph,state->box,state->x);
2344 }
2345
2346 GMX_BARRIER(cr->mpi_comm_mygroup);
2347 GMX_MPE_LOG(ev_update_finish);
2348
2349 if (vsite != NULL)
2350 {
2351 wallcycle_start(wcycle,ewcVSITECONSTR);
2352 if (graph != NULL)
2353 {
2354 shift_self(graph,state->box,state->x);
2355 }
2356 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
2357 top->idef.iparams,top->idef.il,
2358 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
2359
2360 if (graph != NULL)
2361 {
2362 unshift_self(graph,state->box,state->x);
2363 }
2364 wallcycle_stop(wcycle,ewcVSITECONSTR);
2365 }
2366
2367 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
2368 if (ir->nstlist == -1 && bFirstIterate)
2369 {
2370 gs.sig[eglsNABNSB] = nlh.nabnsb;
2371 }
2372 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2373 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2374 constr,
2375 bFirstIterate ? &gs : NULL,(step % gs.nstms == 0),
2376 lastbox,
2377 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2378 cglo_flags
2379 | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
2380 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
2381 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
2382 | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
2383 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
2384 | CGLO_CONSTRAINT
2385 );
2386 if (ir->nstlist == -1 && bFirstIterate)
2387 {
2388 nlh.nabnsb = gs.set[eglsNABNSB];
2389 gs.set[eglsNABNSB] = 0;
2390 }
2391 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
2392 /* ############# END CALC EKIN AND PRESSURE ################# */
2393
2394 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
2395 the virial that should probably be addressed eventually. state->veta has better properies,
2396 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
2397 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
2398
2399 if (bIterations &&
2400 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
2401 trace(shake_vir),&tracevir))
2402 {
2403 break;
2404 }
2405 bFirstIterate = FALSE;
2406 }
2407
2408 update_box(fplog,step,ir,mdatoms,state,graph,f,
2409 ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
2410
2411 /* ################# END UPDATE STEP 2 ################# */
2412 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
2413
2414 /* The coordinates (x) were unshifted in update */
2415 /* if (bFFscan && (shellfc==NULL || bConverged))
2416 {
2417 if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
2418 f,NULL,xcopy,
2419 &(top_global->mols),mdatoms->massT,pres))
2420 {
2421 if (gmx_parallel_env_initialized())
2422 {
2423 gmx_finalize();
2424 }
2425 fprintf(stderr,"\n");
2426 exit(0);
2427 }
2428 }*/
2429 if (!bGStat)
2430 {
2431 /* We will not sum ekinh_old,
2432 * so signal that we still have to do it.
2433 */
2434 bSumEkinhOld = TRUE;
2435 }
2436
2437 /* if (bTCR)
2438 {*/
2439 /* Only do GCT when the relaxation of shells (minimization) has converged,
2440 * otherwise we might be coupling to bogus energies.
2441 * In parallel we must always do this, because the other sims might
2442 * update the FF.
2443 */
2444
2445 /* Since this is called with the new coordinates state->x, I assume
2446 * we want the new box state->box too. / EL 20040121
2447 */
2448 /* do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
2449 ir,MASTER(cr),
2450 mdatoms,&(top->idef),mu_aver,
2451 top_global->mols.nr,cr,
2452 state->box,total_vir,pres,
2453 mu_tot,state->x,f,bConverged);
2454 debug_gmx();
2455 }*/
2456
2457 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
2458
2459 sum_dhdl(enerd,state->lambda,ir);
2460 /* use the directly determined last velocity, not actually the averaged half steps */
2461 if (bTrotter && ir->eI==eiVV)
2462 {
2463 enerd->term[F_EKIN] = last_ekin;
2464 }
2465 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
2466
2467 switch (ir->etc)
2468 {
2469 case etcNO:
2470 break;
2471 case etcBERENDSEN:
2472 break;
2473 case etcNOSEHOOVER:
2474 if (IR_NVT_TROTTER(ir)) {
2475 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + last_conserved;
2476 } else {
2477 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] +
2478 NPT_energy(ir,state,&MassQ);
2479 }
2480 break;
2481 case etcVRESCALE:
2482 enerd->term[F_ECONSERVED] =
2483 enerd->term[F_ETOT] + vrescale_energy(&(ir->opts),
2484 state->therm_integral);
2485 break;
2486 default:
2487 break;
2488 }
2489
2490 /* Check for excessively large energies */
2491 /* if (bIonize)
2492 {
2493 #ifdef GMX_DOUBLE
2494 real etot_max = 1e200;
2495 #else
2496 real etot_max = 1e30;
2497 #endif
2498 if (fabs(enerd->term[F_ETOT]) > etot_max)
2499 {
2500 fprintf(stderr,"Energy too large (%g), giving up\n",
2501 enerd->term[F_ETOT]);
2502 }
2503 }*/
2504 /* ######### END PREPARING EDR OUTPUT ########### */
2505
2506 /* Time for performance */
2507 if (((step % stepout) == 0) || bLastStep)
2508 {
2509 runtime_upd_proc(runtime);
2510 }
2511
2512 /* Output stuff */
2513 if (MASTER(cr))
2514 {
2515 gmx_bool do_dr,do_or;
2516
2517 if (!(bStartingFromCpt && (EI_VV(ir->eI))))
2518 {
2519 if (bNstEner)
2520 {
2521 upd_mdebin(mdebin,bDoDHDL,TRUE,
2522 t,mdatoms->tmass,enerd,state,lastbox,
2523 shake_vir,force_vir,total_vir,pres,
2524 ekind,mu_tot,constr);
2525 }
2526 else
2527 {
2528 upd_mdebin_step(mdebin);
2529 }
2530
2531 do_dr = do_per_step(step,ir->nstdisreout);
2532 do_or = do_per_step(step,ir->nstorireout);
2533
2534 print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
2535 step,t,
2536 eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
2537 }
2538 if (ir->ePull != epullNO)
2539 {
2540 pull_print_output(ir->pull,step,t);
2541 }
2542
2543 if (do_per_step(step,ir->nstlog))
2544 {
2545 if(fflush(fplog) != 0)
2546 {
2547 gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of disk space?");
2548 }
2549 }
2550 }
2551
2552
2553 /* Remaining runtime */
2554 if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
2555 {
2556 if (shellfc)
2557 {
2558 fprintf(stderr,"\n");
2559 }
2560 print_time(stderr,runtime,step,ir,cr);
2561 }
2562
2563 /* Set new positions for the group to embed */
2564 if(!bLastStep){
2565 if(step_rel<=it_xy)
2566 {
2567 fac[0]+=xy_step;
2568 fac[1]+=xy_step;
2569 } else if (step_rel<=(it_xy+it_z))
2570 {
2571 fac[2]+=z_step;
2572 }
2573 resize(ins_at,r_ins,state_global->x,pos_ins,fac);
2574 }
2575
2576 /* Replica exchange */
2577 /* bExchanged = FALSE;
2578 if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
2579 do_per_step(step,repl_ex_nst))
2580 {
2581 bExchanged = replica_exchange(fplog,cr,repl_ex,
2582 state_global,enerd->term,
2583 state,step,t);
2584 }
2585 if (bExchanged && PAR(cr))
2586 {
2587 if (DOMAINDECOMP(cr))
2588 {
2589 dd_partition_system(fplog,step,cr,TRUE,1,
2590 state_global,top_global,ir,
2591 state,&f,mdatoms,top,fr,
2592 vsite,shellfc,constr,
2593 nrnb,wcycle,FALSE);
2594 }
2595 else
2596 {
2597 bcast_state(cr,state,FALSE);
2598 }
2599 }*/
2600
2601 bFirstStep = FALSE;
2602 bInitStep = FALSE;
2603 bStartingFromCpt = FALSE;
2604
2605 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
2606 /* With all integrators, except VV, we need to retain the pressure
2607 * at the current step for coupling at the next step.
2608 */
2609 if ((state->flags & (1<<estPRES_PREV)) &&
2610 (bGStatEveryStep ||
2611 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
2612 {
2613 /* Store the pressure in t_state for pressure coupling
2614 * at the next MD step.
2615 */
2616 copy_mat(pres,state->pres_prev);
2617 }
2618
2619 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
2620
2621 if (bRerunMD)
2622 {
2623 /* read next frame from input trajectory */
2624 bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
2625 }
2626
2627 if (!bRerunMD || !rerun_fr.bStep)
2628 {
2629 /* increase the MD step number */
2630 step++;
2631 step_rel++;
2632 }
2633
2634 cycles = wallcycle_stop(wcycle,ewcSTEP);
2635 if (DOMAINDECOMP(cr) && wcycle)
2636 {
2637 dd_cycles_add(cr->dd,cycles,ddCyclStep);
2638 }
2639
2640 if (step_rel == wcycle_get_reset_counters(wcycle) ||
2641 gs.set[eglsRESETCOUNTERS] != 0)
2642 {
2643 /* Reset all the counters related to performance over the run */
2644 reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime);
2645 wcycle_set_reset_counters(wcycle,-1);
2646 bResetCountersHalfMaxH = FALSE;
2647 gs.set[eglsRESETCOUNTERS] = 0;
2648 }
2649 }
2650 /* End of main MD loop */
2651 debug_gmx();
2652 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
2653 *top_global->name,top_global,
2654 state_global->x,state_global->v,
2655 ir->ePBC,state->box);
2656
2657 /* Stop the time */
2658 runtime_end(runtime);
2659
2660 if (bRerunMD)
2661 {
2662 close_trj(status);
2663 }
2664
2665 if (!(cr->duty & DUTY_PME))
2666 {
2667 /* Tell the PME only node to finish */
2668 gmx_pme_finish(cr);
2669 }
2670
2671 if (MASTER(cr))
2672 {
2673 if (ir->nstcalcenergy > 0 && !bRerunMD)
2674 {
2675 print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
2676 eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
2677 }
2678 }
2679
2680 done_mdoutf(outf);
2681
2682 debug_gmx();
2683
2684 if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
2685 {
2686 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)));
2687 fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
2688 }
2689
2690 if (shellfc && fplog)
2691 {
2692 fprintf(fplog,"Fraction of iterations that converged: %.2f %%\n",
2693 (nconverged*100.0)/step_rel);
2694 fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
2695 tcount/step_rel);
2696 }
2697
2698 /* if (repl_ex_nst > 0 && MASTER(cr))
2699 {
2700 print_replica_exchange_statistics(fplog,repl_ex);
2701 }*/
2702
2703 runtime->nsteps_done = step_rel;
2704
2705 return 0;
2706 }
2707
2708
2709 int mdrunner_membed(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
2710 const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
2711 int nstglobalcomm,
2712 ivec ddxyz,int dd_node_order,real rdd,real rconstr,
2713 const char *dddlb_opt,real dlb_scale,
2714 const char *ddcsx,const char *ddcsy,const char *ddcsz,
2715 int nstepout,int resetstep,int nmultisim,int repl_ex_nst,int repl_ex_seed,
2716 real pforce,real cpt_period,real max_hours,
2717 const char *deviceOptions,
2718 unsigned long Flags,
2719 real xy_fac, real xy_max, real z_fac, real z_max,
2720 int it_xy, int it_z, real probe_rad, int low_up_rm,
2721 int pieces, gmx_bool bALLOW_ASYMMETRY, int maxwarn)
2722 {
2723 double nodetime=0,realtime;
2724 t_inputrec *inputrec;
2725 t_state *state=NULL;
2726 matrix box;
2727 gmx_ddbox_t ddbox;
2728 int npme_major,npme_minor;
2729 real tmpr1,tmpr2;
2730 t_nrnb *nrnb;
2731 gmx_mtop_t *mtop=NULL;
2732 t_mdatoms *mdatoms=NULL;
2733 t_forcerec *fr=NULL;
2734 t_fcdata *fcd=NULL;
2735 real ewaldcoeff=0;
2736 gmx_pme_t *pmedata=NULL;
2737 gmx_vsite_t *vsite=NULL;
2738 gmx_constr_t constr;
2739 int i,m,nChargePerturbed=-1,status,nalloc;
2740 char *gro;
2741 gmx_wallcycle_t wcycle;
2742 gmx_bool bReadRNG,bReadEkin;
2743 int list;
2744 gmx_runtime_t runtime;
2745 int rc;
2746 gmx_large_int_t reset_counters;
2747 gmx_edsam_t ed=NULL;
2748 t_commrec *cr_old=cr;
2749 int nthreads=1,nthreads_requested=1;
2750
2751
2752 char *ins;
2753 int rm_bonded_at,fr_id,fr_i=0,tmp_id,warn=0;
2754 int ng,j,max_lip_rm,ins_grp_id,ins_nat,mem_nat,ntype,lip_rm,tpr_version;
2755 real xy_step=0,z_step=0;
2756 real prot_area;
2757 rvec *r_ins=NULL,fac;
2758 t_block *ins_at,*rest_at;
2759 pos_ins_t *pos_ins;
2760 mem_t *mem_p;
2761 rmm_t *rm_p;
2762 gmx_groups_t *groups;
2763 gmx_bool bExcl=FALSE;
2764 t_atoms atoms;
2765 t_pbc *pbc;
2766 char **piecename=NULL;
2767
2768 /* CAUTION: threads may be started later on in this function, so
2769 cr doesn't reflect the final parallel state right now */
2770 snew(inputrec,1);
2771 snew(mtop,1);
2772
2773 if (bVerbose && SIMMASTER(cr))
2774 {
2775 fprintf(stderr,"Getting Loaded...\n");
2776 }
2777
2778 if (Flags & MD_APPENDFILES)
2779 {
2780 fplog = NULL;
2781 }
2782
2783 snew(state,1);
2784 if (MASTER(cr))
2785 {
2786 /* Read (nearly) all data required for the simulation */
2787 read_tpx_state(ftp2fn(efTPX,nfile,fnm),inputrec,state,NULL,mtop);
2788
2789 /* NOW the threads will be started: */
2790 #ifdef GMX_THREADS
2791 #endif
2792 }
2793 /* END OF CAUTION: cr is now reliable */
2794
2795 if (PAR(cr))
2796 {
2797 /* now broadcast everything to the non-master nodes/threads: */
2798 init_parallel(fplog, cr, inputrec, mtop);
2799 }
2800 /* now make sure the state is initialized and propagated */
2801 set_state_entries(state,inputrec,cr->nnodes);
2802
2803 if (can_use_allvsall(inputrec,mtop,TRUE,cr,fplog))
2804 {
2805 /* All-vs-all loops do not work with domain decomposition */
2806 Flags |= MD_PARTDEC;
2807 }
2808
2809 if (!EEL_PME(inputrec->coulombtype) || (Flags & MD_PARTDEC))
2810 {
2811 cr->npmenodes = 0;
2812 }
2813
2814 snew(ins_at,1);
2815 snew(pos_ins,1);
2816 if(MASTER(cr))
2817 {
2818 tpr_version = get_tpr_version(ftp2fn(efTPX,nfile,fnm));
2819 if (tpr_version<58)
2820 gmx_fatal(FARGS,"Version of *.tpr file to old (%d). Rerun grompp with gromacs VERSION 4.0.3 or newer.\n",tpr_version);
2821
2822 if( inputrec->eI != eiMD )
2823 gmx_input("Change integrator to md in mdp file.");
2824
2825 if(PAR(cr))
2826 gmx_input("Sorry, parallel g_membed is not yet fully functrional.");
2827
2828 groups=&(mtop->groups);
2829
2830 atoms=gmx_mtop_global_atoms(mtop);
2831 snew(mem_p,1);
2832 fprintf(stderr,"\nSelect a group to embed in the membrane:\n");
2833 get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),1,&(ins_at->nr),&(ins_at->index),&ins);
2834 ins_grp_id = search_string(ins,groups->ngrpname,(groups->grpname));
2835 fprintf(stderr,"\nSelect a group to embed %s into (e.g. the membrane):\n",ins);
2836 get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),1,&(mem_p->mem_at.nr),&(mem_p->mem_at.index),&(mem_p->name));
2837
2838 pos_ins->pieces=pieces;
2839 snew(pos_ins->nidx,pieces);
2840 snew(pos_ins->subindex,pieces);
2841 snew(piecename,pieces);
2842 if (pieces>1)
2843 {
2844 fprintf(stderr,"\nSelect pieces to embed:\n");
2845 get_index(&atoms,ftp2fn_null(efNDX,nfile,fnm),pieces,pos_ins->nidx,pos_ins->subindex,piecename);
2846 }
2847 else
2848 {
2849 /*use whole embedded group*/
2850 snew(pos_ins->nidx,1);
2851 snew(pos_ins->subindex,1);
2852 pos_ins->nidx[0]=ins_at->nr;
2853 pos_ins->subindex[0]=ins_at->index;
2854 }
2855
2856 if(probe_rad<0.2199999)
2857 {
2858 warn++;
2859 fprintf(stderr,"\nWarning %d:\nA probe radius (-rad) smaller than 0.2 can result in overlap between waters "
2860 "and the group to embed, which will result in Lincs errors etc.\nIf you are sure, you can increase maxwarn.\n\n",warn);
2861 }
2862
2863 if(xy_fac<0.09999999)
2864 {
2865 warn++;
2866 fprintf(stderr,"\nWarning %d:\nThe initial size of %s is probably too smal.\n"
2867 "If you are sure, you can increase maxwarn.\n\n",warn,ins);
2868 }
2869
2870 if(it_xy<1000)
2871 {
2872 warn++;
2873 fprintf(stderr,"\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d) is probably too small.\n"
2874 "Increase -nxy or, if you are sure, you can increase maxwarn.\n\n",warn,ins,it_xy);
2875 }
2876
2877 if( (it_z<100) && ( z_fac<0.99999999 || z_fac>1.0000001) )
2878 {
2879 warn++;
2880 fprintf(stderr,"\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d) is probably too small.\n"
2881 "Increase -nz or, if you are sure, you can increase maxwarn.\n\n",warn,ins,it_z);
2882 }
2883
2884 if(it_xy+it_z>inputrec->nsteps)
2885 {
2886 warn++;
2887 fprintf(stderr,"\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the number of steps in the tpr.\n"
2888 "If you are sure, you can increase maxwarn.\n\n",warn);
2889 }
2890
2891 fr_id=-1;
2892 if( inputrec->opts.ngfrz==1)
2893 gmx_fatal(FARGS,"You did not specify \"%s\" as a freezegroup.",ins);
2894 for(i=0;i<inputrec->opts.ngfrz;i++)
2895 {
2896 tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i];
2897 if(ins_grp_id==tmp_id)
2898 {
2899 fr_id=tmp_id;
2900 fr_i=i;
2901 }
2902 }
2903 if (fr_id == -1 )
2904 gmx_fatal(FARGS,"\"%s\" not as freezegroup defined in the mdp-file.",ins);
2905
2906 for(i=0;i<DIM;i++)
2907 if( inputrec->opts.nFreeze[fr_i][i] != 1)
2908 gmx_fatal(FARGS,"freeze dimensions for %s are not Y Y Y\n",ins);
2909
2910 ng = groups->grps[egcENER].nr;
2911 if (ng == 1)
2912 gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
2913
2914 for(i=0;i<ng;i++)
2915 {
2916 for(j=0;j<ng;j++)
2917 {
2918 if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL)
2919 {
2920 bExcl = TRUE;
2921 if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) || (groups->grps[egcENER].nm_ind[j] != ins_grp_id) )
2922 gmx_fatal(FARGS,"Energy exclusions \"%s\" and \"%s\" do not match the group to embed \"%s\"",
2923 *groups->grpname[groups->grps[egcENER].nm_ind[i]],
2924 *groups->grpname[groups->grps[egcENER].nm_ind[j]],ins);
2925 }
2926 }
2927 }
2928 if (!bExcl)
2929 gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in the freeze group");
2930
2931 /* Set all atoms in box*/
2932 /*set_inbox(state->natoms,state->x);*/
2933
2934 /* Guess the area the protein will occupy in the membrane plane Calculate area per lipid*/
2935 snew(rest_at,1);
2936 ins_nat = init_ins_at(ins_at,rest_at,state,pos_ins,groups,ins_grp_id,xy_max);
2937 /* Check moleculetypes in insertion group */
2938 check_types(ins_at,rest_at,mtop);
2939
2940 mem_nat = init_mem_at(mem_p,mtop,state->x,state->box,pos_ins);
2941
2942 prot_area = est_prot_area(pos_ins,state->x,ins_at,mem_p);
2943 if ( (prot_area>7.5) && ( (state->box[XX][XX]*state->box[YY][YY]-state->box[XX][YY]*state->box[YY][XX])<50) )
2944 {
2945 warn++;
2946 fprintf(stderr,"\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
2947 "This might cause pressure problems during the growth phase. Just try with\n"
2948 "current setup (-maxwarn + 1), but if pressure problems occur, lower the\n"
2949 "compressibility in the mdp-file or use no pressure coupling at all.\n\n",warn);
2950 }
2951 if(warn>maxwarn)
2952 gmx_fatal(FARGS,"Too many warnings.\n");
2953
2954 printf("The estimated area of the protein in the membrane is %.3f nm^2\n",prot_area);
2955 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);
2956
2957 /* Maximum number of lipids to be removed*/
2958 max_lip_rm=(int)(2*prot_area/mem_p->lip_area);
2959 printf("Maximum number of lipids that will be removed is %d.\n",max_lip_rm);
2960
2961 printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
2962 "This resizing will be done with respect to the geometrical center of all protein atoms\n"
2963 "that span the membrane region, i.e. z between %.3f and %.3f\n\n",xy_fac,z_fac,mem_p->zmin,mem_p->zmax);
2964
2965 /* resize the protein by xy and by z if necessary*/
2966 snew(r_ins,ins_at->nr);
2967 init_resize(ins_at,r_ins,pos_ins,mem_p,state->x,bALLOW_ASYMMETRY);
2968 fac[0]=fac[1]=xy_fac;
2969 fac[2]=z_fac;
2970
2971 xy_step =(xy_max-xy_fac)/(double)(it_xy);
2972 z_step =(z_max-z_fac)/(double)(it_z-1);
2973
2974 resize(ins_at,r_ins,state->x,pos_ins,fac);
2975
2976 /* remove overlapping lipids and water from the membrane box*/
2977 /*mark molecules to be removed*/
2978 snew(pbc,1);
2979 set_pbc(pbc,inputrec->ePBC,state->box);
2980
2981 snew(rm_p,1);
2982 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);
2983 lip_rm -= low_up_rm;
2984
2985 if(fplog)
2986 for(i=0;i<rm_p->nr;i++)
2987 fprintf(fplog,"rm mol %d\n",rm_p->mol[i]);
2988
2989 for(i=0;i<mtop->nmolblock;i++)
2990 {
2991 ntype=0;
2992 for(j=0;j<rm_p->nr;j++)
2993 if(rm_p->block[j]==i)
2994 ntype++;
2995 printf("Will remove %d %s molecules\n",ntype,*(mtop->moltype[mtop->molblock[i].type].name));
2996 }
2997
2998 if(lip_rm>max_lip_rm)
2999 {
3000 warn++;
3001 fprintf(stderr,"\nWarning %d:\nTrying to remove a larger lipid area than the estimated protein area\n"
3002 "Try making the -xyinit resize factor smaller. If you are sure about this increase maxwarn.\n\n",warn);
3003 }
3004
3005 /*remove all lipids and waters overlapping and update all important structures*/
3006 rm_group(inputrec,groups,mtop,rm_p,state,ins_at,pos_ins);
3007
3008 rm_bonded_at = rm_bonded(ins_at,mtop);
3009 if (rm_bonded_at != ins_at->nr)
3010 {
3011 fprintf(stderr,"Warning: The number of atoms for which the bonded interactions are removed is %d, "
3012 "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
3013 "molecule type as atoms that are not to be embedded.\n",rm_bonded_at,ins_at->nr);
3014 }
3015
3016 if(warn>maxwarn)
3017 gmx_fatal(FARGS,"Too many warnings.\nIf you are sure these warnings are harmless, you can increase -maxwarn");
3018
3019 if (MASTER(cr))
3020 {
3021 if (ftp2bSet(efTOP,nfile,fnm))
3022 top_update(opt2fn("-p",nfile,fnm),ins,rm_p,mtop);
3023 }
3024
3025 sfree(pbc);
3026 sfree(rest_at);
3027 }
3028
3029 #ifdef GMX_FAHCORE
3030 fcRegisterSteps(inputrec->nsteps,inputrec->init_step);
3031 #endif
3032
3033 /* NMR restraints must be initialized before load_checkpoint,
3034 * since with time averaging the history is added to t_state.
3035 * For proper consistency check we therefore need to extend
3036 * t_state here.
3037 * So the PME-only nodes (if present) will also initialize
3038 * the distance restraints.
3039 */
3040 snew(fcd,1);
3041
3042 /* This needs to be called before read_checkpoint to extend the state */
3043 init_disres(fplog,mtop,inputrec,cr,Flags & MD_PARTDEC,fcd,state,FALSE);
3044
3045 if (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0)
3046 {
3047 if (PAR(cr) && !(Flags & MD_PARTDEC))
3048 {
3049 gmx_fatal(FARGS,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
3050 }
3051 /* Orientation restraints */
3052 if (MASTER(cr))
3053 {
3054 init_orires(fplog,mtop,state->x,inputrec,cr->ms,&(fcd->orires),
3055 state);
3056 }
3057 }
3058
3059 if (DEFORM(*inputrec))
3060 {
3061 /* Store the deform reference box before reading the checkpoint */
3062 if (SIMMASTER(cr))
3063 {
3064 copy_mat(state->box,box);
3065 }
3066 if (PAR(cr))
3067 {
3068 gmx_bcast(sizeof(box),box,cr);
3069 }
3070 /* Because we do not have the update struct available yet
3071 * in which the reference values should be stored,
3072 * we store them temporarily in static variables.
3073 * This should be thread safe, since they are only written once
3074 * and with identical values.
3075 */
3076 /* deform_init_init_step_tpx = inputrec->init_step;*/
3077 /* copy_mat(box,deform_init_box_tpx);*/
3078 }
3079
3080 if (opt2bSet("-cpi",nfile,fnm))
3081 {
3082 /* Check if checkpoint file exists before doing continuation.
3083 * This way we can use identical input options for the first and subsequent runs...
3084 */
3085 if( gmx_fexist_master(opt2fn_master("-cpi",nfile,fnm,cr),cr) )
3086 {
3087 load_checkpoint(opt2fn_master("-cpi",nfile,fnm,cr),&fplog,
3088 cr,Flags & MD_PARTDEC,ddxyz,
3089 inputrec,state,&bReadRNG,&bReadEkin,
3090 (Flags & MD_APPENDFILES),
3091 (Flags & MD_APPENDFILESSET));
3092
3093 if (bReadRNG)
3094 {
3095 Flags |= MD_READ_RNG;
3096 }
3097 if (bReadEkin)
3098 {
3099 Flags |= MD_READ_EKIN;
3100 }
3101 }
3102 }
3103
3104 if ((MASTER(cr) || (Flags & MD_SEPPOT)) && (Flags & MD_APPENDFILES))
3105 {
3106 gmx_log_open(ftp2fn(efLOG,nfile,fnm),cr,!(Flags & MD_SEPPOT),
3107 Flags,&fplog);
3108 }
3109
3110 if (SIMMASTER(cr))
3111 {
3112 copy_mat(state->box,box);
3113 }
3114
3115 if (PAR(cr))
3116 {
3117 gmx_bcast(sizeof(box),box,cr);
3118 }
3119
3120 if (bVerbose && SIMMASTER(cr))
3121 {
3122 fprintf(stderr,"Loaded with Money\n\n");
3123 }
3124
3125 if (PAR(cr) && !((Flags & MD_PARTDEC) || EI_TPI(inputrec->eI)))
3126 {
3127 cr->dd = init_domain_decomposition(fplog,cr,Flags,ddxyz,rdd,rconstr,
3128 dddlb_opt,dlb_scale,
3129 ddcsx,ddcsy,ddcsz,
3130 mtop,inputrec,
3131 box,state->x,
3132 &ddbox,&npme_major,&npme_minor);
3133
3134 make_dd_communicators(fplog,cr,dd_node_order);
3135
3136 /* Set overallocation to avoid frequent reallocation of arrays */
3137 set_over_alloc_dd(TRUE);
3138 }
3139 else
3140 {
3141 /* PME, if used, is done on all nodes with 1D decomposition */
3142 cr->npmenodes = 0;
3143 cr->duty = (DUTY_PP | DUTY_PME);
3144 npme_major = cr->nnodes;
3145 npme_minor = 1;
3146
3147 if (inputrec->ePBC == epbcSCREW)
3148 {
3149 gmx_fatal(FARGS,
3150 "pbc=%s is only implemented with domain decomposition",
3151 epbc_names[inputrec->ePBC]);
3152 }
3153 }
3154
3155 if (PAR(cr))
3156 {
3157 /* After possible communicator splitting in make_dd_communicators.
3158 * we can set up the intra/inter node communication.
3159 */
3160 gmx_setup_nodecomm(fplog,cr);
3161 }
3162
3163 wcycle = wallcycle_init(fplog,resetstep,cr);
3164 if (PAR(cr))
3165 {
3166 /* Master synchronizes its value of reset_counters with all nodes
3167 * including PME only nodes */
3168 reset_counters = wcycle_get_reset_counters(wcycle);
3169 gmx_bcast_sim(sizeof(reset_counters),&reset_counters,cr);
3170 wcycle_set_reset_counters(wcycle, reset_counters);
3171 }
3172
3173
3174 snew(nrnb,1);
3175 if (cr->duty & DUTY_PP)
3176 {
3177 /* For domain decomposition we allocate dynamically
3178 * in dd_partition_system.
3179 */
3180 if (DOMAINDECOMP(cr))
3181 {
3182 bcast_state_setup(cr,state);
3183 }
3184 else
3185 {
3186 if (PAR(cr))
3187 {
3188 if (!MASTER(cr))
3189 {
3190 snew(state,1);
3191 }
3192 bcast_state(cr,state,TRUE);
3193 }
3194 }
3195
3196 /* Dihedral Restraints */
3197 if (gmx_mtop_ftype_count(mtop,F_DIHRES) > 0)
3198 {
3199 init_dihres(fplog,mtop,inputrec,fcd);
3200 }
3201
3202 /* Initiate forcerecord */
3203 fr = mk_forcerec();
3204 init_forcerec(fplog,oenv,fr,fcd,inputrec,mtop,cr,box,FALSE,
3205 opt2fn("-table",nfile,fnm),
3206 opt2fn("-tablep",nfile,fnm),
3207 opt2fn("-tableb",nfile,fnm),FALSE,pforce);
3208
3209 /* version for PCA_NOT_READ_NODE (see md.c) */
3210 /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
3211 "nofile","nofile","nofile",FALSE,pforce);
3212 */
3213 fr->bSepDVDL = ((Flags & MD_SEPPOT) == MD_SEPPOT);
3214
3215 /* Initialize QM-MM */
3216 if(fr->bQMMM)
3217 {
3218 init_QMMMrec(cr,box,mtop,inputrec,fr);
3219 }
3220
3221 /* Initialize the mdatoms structure.
3222 * mdatoms is not filled with atom data,
3223 * as this can not be done now with domain decomposition.
3224 */
3225 mdatoms = init_mdatoms(fplog,mtop,inputrec->efep!=efepNO);
3226
3227 /* Initialize the virtual site communication */
3228 vsite = init_vsite(mtop,cr);
3229
3230 calc_shifts(box,fr->shift_vec);
3231
3232 /* With periodic molecules the charge groups should be whole at start up
3233 * and the virtual sites should not be far from their proper positions.
3234 */
3235 if (!inputrec->bContinuation && MASTER(cr) &&
3236 !(inputrec->ePBC != epbcNONE && inputrec->bPeriodicMols))
3237 {
3238 /* Make molecules whole at start of run */
3239 if (fr->ePBC != epbcNONE)
3240 {
3241 do_pbc_first_mtop(fplog,inputrec->ePBC,box,mtop,state->x);
3242 }
3243 if (vsite)
3244 {
3245 /* Correct initial vsite positions are required
3246 * for the initial distribution in the domain decomposition
3247 * and for the initial shell prediction.
3248 */
3249 construct_vsites_mtop(fplog,vsite,mtop,state->x);
3250 }
3251 }
3252
3253 /* Initiate PPPM if necessary */
3254 if (fr->eeltype == eelPPPM)
3255 {
3256 if (mdatoms->nChargePerturbed)
3257 {
3258 gmx_fatal(FARGS,"Free energy with %s is not implemented",
3259 eel_names[fr->eeltype]);
3260 }
3261 status = gmx_pppm_init(fplog,cr,oenv,FALSE,TRUE,box,
3262 getenv("GMXGHAT"),inputrec, (Flags & MD_REPRODUCIBLE));
3263 if (status != 0)
3264 {
3265 gmx_fatal(FARGS,"Error %d initializing PPPM",status);
3266 }
3267 }
3268
3269 if (EEL_PME(fr->eeltype))
3270 {
3271 ewaldcoeff = fr->ewaldcoeff;
3272 pmedata = &fr->pmedata;
3273 }
3274 else
3275 {
3276 pmedata = NULL;
3277 }
3278 }
3279 else
3280 {
3281 /* This is a PME only node */
3282
3283 /* We don't need the state */
3284 done_state(state);
3285
3286 ewaldcoeff = calc_ewaldcoeff(inputrec->rcoulomb, inputrec->ewald_rtol);
3287 snew(pmedata,1);
3288 }
3289
3290 /* Initiate PME if necessary,
3291 * either on all nodes or on dedicated PME nodes only. */
3292 if (EEL_PME(inputrec->coulombtype))
3293 {
3294 if (mdatoms)
3295 {
3296 nChargePerturbed = mdatoms->nChargePerturbed;
3297 }
3298 if (cr->npmenodes > 0)
3299 {
3300 /* The PME only nodes need to know nChargePerturbed */
3301 gmx_bcast_sim(sizeof(nChargePerturbed),&nChargePerturbed,cr);
3302 }
3303 if (cr->duty & DUTY_PME)
3304 {
3305 status = gmx_pme_init(pmedata,cr,npme_major,npme_minor,inputrec,
3306 mtop ? mtop->natoms : 0,nChargePerturbed,
3307 (Flags & MD_REPRODUCIBLE));
3308 if (status != 0)
3309 {
3310 gmx_fatal(FARGS,"Error %d initializing PME",status);
3311 }
3312 }
3313 }
3314
3315
3316 /* if (integrator[inputrec->eI].func == do_md
3317 #ifdef GMX_OPENMM
3318 ||
3319 integrator[inputrec->eI].func == do_md_openmm
3320 #endif
3321 )
3322 {*/
3323 /* Turn on signal handling on all nodes */
3324 /*
3325 * (A user signal from the PME nodes (if any)
3326 * is communicated to the PP nodes.
3327 */
3328 signal_handler_install();
3329 /* }*/
3330
3331 if (cr->duty & DUTY_PP)
3332 {
3333 if (inputrec->ePull != epullNO)
3334 {
3335 /* Initialize pull code */
3336 init_pull(fplog,inputrec,nfile,fnm,mtop,cr,oenv,
3337 EI_DYNAMICS(inputrec->eI) && MASTER(cr),Flags);
3338 }
3339
3340 constr = init_constraints(fplog,mtop,inputrec,ed,state,cr);
3341
3342 if (DOMAINDECOMP(cr))
3343 {
3344 dd_init_bondeds(fplog,cr->dd,mtop,vsite,constr,inputrec,
3345 Flags & MD_DDBONDCHECK,fr->cginfo_mb);
3346
3347 set_dd_parameters(fplog,cr->dd,dlb_scale,inputrec,fr,&ddbox);
3348
3349 setup_dd_grid(fplog,cr->dd);
3350 }
3351
3352 /* Now do whatever the user wants us to do (how flexible...) */
3353 do_md_membed(fplog,cr,nfile,fnm,
3354 oenv,bVerbose,bCompact,
3355 nstglobalcomm,
3356 vsite,constr,
3357 nstepout,inputrec,mtop,
3358 fcd,state,
3359 mdatoms,nrnb,wcycle,ed,fr,
3360 repl_ex_nst,repl_ex_seed,
3361 cpt_period,max_hours,
3362 deviceOptions,
3363 Flags,
3364 &runtime,
3365 fac, r_ins, pos_ins, ins_at,
3366 xy_step, z_step, it_xy, it_z);
3367
3368 if (inputrec->ePull != epullNO)
3369 {
3370 finish_pull(fplog,inputrec->pull);
3371 }
3372 }
3373 else
3374 {
3375 /* do PME only */
3376 gmx_pmeonly(*pmedata,cr,nrnb,wcycle,ewaldcoeff,FALSE,inputrec);
3377 }
3378
3379 if (EI_DYNAMICS(inputrec->eI) || EI_TPI(inputrec->eI))
3380 {
3381 /* Some timing stats */
3382 if (MASTER(cr))
3383 {
3384 if (runtime.proc == 0)
3385 {
3386 runtime.proc = runtime.real;
3387 }
3388 }
3389 else
3390 {
3391 runtime.real = 0;
3392 }
3393 }
3394
3395 wallcycle_stop(wcycle,ewcRUN);
3396
3397 /* Finish up, write some stuff
3398 * if rerunMD, don't write last frame again
3399 */
3400 finish_run(fplog,cr,ftp2fn(efSTO,nfile,fnm),
3401 inputrec,nrnb,wcycle,&runtime,
3402 EI_DYNAMICS(inputrec->eI) && !MULTISIM(cr));
3403
3404 /* Does what it says */
3405 print_date_and_time(fplog,cr->nodeid,"Finished mdrun",&runtime);
3406
3407 /* Close logfile already here if we were appending to it */
3408 if (MASTER(cr) && (Flags & MD_APPENDFILES))
3409 {
3410 gmx_log_close(fplog);
3411 }
3412
3413 if (pieces>1)
3414 {
3415 sfree(piecename);
3416 }
3417
3418 rc=(int)gmx_get_stop_condition();
3419
3420 return rc;
3421 }
3422
3423 int gmx_membed(int argc,char *argv[])
3424 {
3425 const char *desc[] = {
3426 "[TT]g_membed[tt] embeds a membrane protein into an equilibrated lipid bilayer at the position",
3427 "and orientation specified by the user.[PAR]",
3428 "SHORT MANUAL[BR]------------[BR]",
3429 "The user should merge the structure files of the protein and membrane (+solvent), creating a",
3430 "single structure file with the protein overlapping the membrane at the desired position and",
3431 "orientation. The box size is taken from the membrane structure file. The corresponding topology",
3432 "files should also be merged. Consecutively, create a [TT].tpr[tt] file (input for [TT]g_membed[tt]) from these files,"
3433 "with the following options included in the [TT].mdp[tt] file.[BR]",
3434 " - [TT]integrator = md[tt][BR]",
3435 " - [TT]energygrp = Protein[tt] (or other group that you want to insert)[BR]",
3436 " - [TT]freezegrps = Protein[tt][BR]",
3437 " - [TT]freezedim = Y Y Y[tt][BR]",
3438 " - [TT]energygrp_excl = Protein Protein[tt][BR]",
3439 "The output is a structure file containing the protein embedded in the membrane. If a topology",
3440 "file is provided, the number of lipid and ",
3441 "solvent molecules will be updated to match the new structure file.[BR]",
3442 "For a more extensive manual see Wolf et al, J Comp Chem 31 (2010) 2169-2174, Appendix.[PAR]",
3443 "SHORT METHOD DESCRIPTION[BR]",
3444 "------------------------[BR]",
3445 "1. The protein is resized around its center of mass by a factor [TT]-xy[tt] in the xy-plane",
3446 "(the membrane plane) and a factor [TT]-z[tt] in the [IT]z[it]-direction (if the size of the",
3447 "protein in the z-direction is the same or smaller than the width of the membrane, a",
3448 "[TT]-z[tt] value larger than 1 can prevent that the protein will be enveloped by the lipids).[BR]",
3449 "2. All lipid and solvent molecules overlapping with the resized protein are removed. All",
3450 "intraprotein interactions are turned off to prevent numerical issues for small values of [TT]-xy[tt]",
3451 " or [TT]-z[tt][BR]",
3452 "3. One md step is performed.[BR]",
3453 "4. The resize factor ([TT]-xy[tt] or [TT]-z[tt]) is incremented by a small amount ((1-xy)/nxy or (1-z)/nz) and the",
3454 "protein is resized again around its center of mass. The resize factor for the xy-plane",
3455 "is incremented first. The resize factor for the z-direction is not changed until the [TT]-xy[tt] factor",
3456 "is 1 (thus after [TT]-nxy[tt] iterations).[BR]",
3457 "5. Repeat step 3 and 4 until the protein reaches its original size ([TT]-nxy[tt] + [TT]-nz[tt] iterations).[BR]",
3458 "For a more extensive method description see Wolf et al, J Comp Chem, 31 (2010) 2169-2174.[PAR]",
3459 "NOTE[BR]----[BR]",
3460 " - Protein can be any molecule you want to insert in the membrane.[BR]",
3461 " - It is recommended to perform a short equilibration run after the embedding",
3462 "(see Wolf et al, J Comp Chem 31 (2010) 2169-2174), to re-equilibrate the membrane. Clearly",
3463 "protein equilibration might require longer.[PAR]"
3464 };
3465 t_commrec *cr;
3466 t_filenm fnm[] = {
3467 { efTPX, "-f", "into_mem", ffREAD },
3468 { efNDX, "-n", "index", ffOPTRD },
3469 { efTOP, "-p", "topol", ffOPTRW },
3470 { efTRN, "-o", NULL, ffWRITE },
3471 { efXTC, "-x", NULL, ffOPTWR },
3472 { efCPT, "-cpi", NULL, ffOPTRD },
3473 { efCPT, "-cpo", NULL, ffOPTWR },
3474 { efSTO, "-c", "membedded", ffWRITE },
3475 { efEDR, "-e", "ener", ffWRITE },
3476 { efLOG, "-g", "md", ffWRITE },
3477 { efEDI, "-ei", "sam", ffOPTRD },
3478 { efTRX, "-rerun", "rerun", ffOPTRD },
3479 { efXVG, "-table", "table", ffOPTRD },
3480 { efXVG, "-tablep", "tablep", ffOPTRD },
3481 { efXVG, "-tableb", "table", ffOPTRD },
3482 { efXVG, "-dhdl", "dhdl", ffOPTWR },
3483 { efXVG, "-field", "field", ffOPTWR },
3484 { efXVG, "-table", "table", ffOPTRD },
3485 { efXVG, "-tablep", "tablep", ffOPTRD },
3486 { efXVG, "-tableb", "table", ffOPTRD },
3487 { efTRX, "-rerun", "rerun", ffOPTRD },
3488 { efXVG, "-tpi", "tpi", ffOPTWR },
3489 { efXVG, "-tpid", "tpidist", ffOPTWR },
3490 { efEDI, "-ei", "sam", ffOPTRD },
3491 { efEDO, "-eo", "sam", ffOPTWR },
3492 { efGCT, "-j", "wham", ffOPTRD },
3493 { efGCT, "-jo", "bam", ffOPTWR },
3494 { efXVG, "-ffout", "gct", ffOPTWR },
3495 { efXVG, "-devout", "deviatie", ffOPTWR },
3496 { efXVG, "-runav", "runaver", ffOPTWR },
3497 { efXVG, "-px", "pullx", ffOPTWR },
3498 { efXVG, "-pf", "pullf", ffOPTWR },
3499 { efMTX, "-mtx", "nm", ffOPTWR },
3500 { efNDX, "-dn", "dipole", ffOPTWR },
3501 { efRND, "-multidir",NULL, ffOPTRDMULT}
3502 };
3503 #define NFILE asize(fnm)
3504
3505 /* Command line options ! */
3506 gmx_bool bCart = FALSE;
3507 gmx_bool bPPPME = FALSE;
3508 gmx_bool bPartDec = FALSE;
3509 gmx_bool bDDBondCheck = TRUE;
3510 gmx_bool bDDBondComm = TRUE;
3511 gmx_bool bVerbose = FALSE;
3512 gmx_bool bCompact = TRUE;
3513 gmx_bool bSepPot = FALSE;
3514 gmx_bool bRerunVSite = FALSE;
3515 gmx_bool bIonize = FALSE;
3516 gmx_bool bConfout = TRUE;
3517 gmx_bool bReproducible = FALSE;
3518
3519 int npme=-1;
3520 int nmultisim=0;
3521 int nstglobalcomm=-1;
3522 int repl_ex_nst=0;
3523 int repl_ex_seed=-1;
3524 int nstepout=100;
3525 int nthreads=0; /* set to determine # of threads automatically */
3526 int resetstep=-1;
3527
3528 rvec realddxyz={0,0,0};
3529 const char *ddno_opt[ddnoNR+1] =
3530 { NULL, "interleave", "pp_pme", "cartesian", NULL };
3531 const char *dddlb_opt[] =
3532 { NULL, "auto", "no", "yes", NULL };
3533 real rdd=0.0,rconstr=0.0,dlb_scale=0.8,pforce=-1;
3534 char *ddcsx=NULL,*ddcsy=NULL,*ddcsz=NULL;
3535 real cpt_period=15.0,max_hours=-1;
3536 gmx_bool bAppendFiles=TRUE,bAddPart=TRUE;
3537 gmx_bool bResetCountersHalfWay=FALSE;
3538 output_env_t oenv=NULL;
3539 const char *deviceOptions = "";
3540
3541 real xy_fac = 0.5;
3542 real xy_max = 1.0;
3543 real z_fac = 1.0;
3544 real z_max = 1.0;
3545 int it_xy = 1000;
3546 int it_z = 0;
3547 real probe_rad = 0.22;
3548 int low_up_rm = 0;
3549 int maxwarn=0;
3550 int pieces=1;
3551 gmx_bool bALLOW_ASYMMETRY=FALSE;
3552
3553
3554 /* arguments relevant to OPENMM only*/
3555 #ifdef GMX_OPENMM
3556 gmx_input("g_membed not functional in openmm");
3557 #endif
3558
3559 t_pargs pa[] = {
3560 { "-xyinit", FALSE, etREAL, {&xy_fac}, "Resize factor for the protein in the xy dimension before starting embedding" },
3561 { "-xyend", FALSE, etREAL, {&xy_max}, "Final resize factor in the xy dimension" },
3562 { "-zinit", FALSE, etREAL, {&z_fac}, "Resize factor for the protein in the z dimension before starting embedding" },
3563 { "-zend", FALSE, etREAL, {&z_max}, "Final resize faction in the z dimension" },
3564 { "-nxy", FALSE, etINT, {&it_xy}, "Number of iteration for the xy dimension" },
3565 { "-nz", FALSE, etINT, {&it_z}, "Number of iterations for the z dimension" },
3566 { "-rad", FALSE, etREAL, {&probe_rad}, "Probe radius to check for overlap between the group to embed and the membrane"},
3567 { "-pieces", FALSE, etINT, {&pieces}, "Perform piecewise resize. Select parts of the group to insert and resize these with respect to their own geometrical center." },
3568 { "-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." },
3569 { "-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." },
3570 { "-maxwarn", FALSE, etINT, {&maxwarn}, "Maximum number of warning allowed" },
3571 { "-pd", FALSE, etBOOL,{&bPartDec},
3572 "HIDDENUse particle decompostion" },
3573 { "-dd", FALSE, etRVEC,{&realddxyz},
3574 "HIDDENDomain decomposition grid, 0 is optimize" },
3575 { "-nt", FALSE, etINT, {&nthreads},
3576 "HIDDENNumber of threads to start (0 is guess)" },
3577 { "-npme", FALSE, etINT, {&npme},
3578 "HIDDENNumber of separate nodes to be used for PME, -1 is guess" },
3579 { "-ddorder", FALSE, etENUM, {ddno_opt},
3580 "HIDDENDD node order" },
3581 { "-ddcheck", FALSE, etBOOL, {&bDDBondCheck},
3582 "HIDDENCheck for all bonded interactions with DD" },
3583 { "-ddbondcomm", FALSE, etBOOL, {&bDDBondComm},
3584 "HIDDENUse special bonded atom communication when [TT]-rdd[tt] > cut-off" },
3585 { "-rdd", FALSE, etREAL, {&rdd},
3586 "HIDDENThe maximum distance for bonded interactions with DD (nm), 0 is determine from initial coordinates" },
3587 { "-rcon", FALSE, etREAL, {&rconstr},
3588 "HIDDENMaximum distance for P-LINCS (nm), 0 is estimate" },
3589 { "-dlb", FALSE, etENUM, {dddlb_opt},
3590 "HIDDENDynamic load balancing (with DD)" },
3591 { "-dds", FALSE, etREAL, {&dlb_scale},
3592 "HIDDENMinimum allowed dlb scaling of the DD cell size" },
3593 { "-ddcsx", FALSE, etSTR, {&ddcsx},
3594 "HIDDENThe DD cell sizes in x" },
3595 { "-ddcsy", FALSE, etSTR, {&ddcsy},
3596 "HIDDENThe DD cell sizes in y" },
3597 { "-ddcsz", FALSE, etSTR, {&ddcsz},
3598 "HIDDENThe DD cell sizes in z" },
3599 { "-gcom", FALSE, etINT,{&nstglobalcomm},
3600 "HIDDENGlobal communication frequency" },
3601 { "-compact", FALSE, etBOOL,{&bCompact},
3602 "Write a compact log file" },
3603 { "-seppot", FALSE, etBOOL, {&bSepPot},
3604 "HIDDENWrite separate V and dVdl terms for each interaction type and node to the log file(s)" },
3605 { "-pforce", FALSE, etREAL, {&pforce},
3606 "HIDDENPrint all forces larger than this (kJ/mol nm)" },
3607 { "-reprod", FALSE, etBOOL,{&bReproducible},
3608 "HIDDENTry to avoid optimizations that affect binary reproducibility" },
3609 { "-multi", FALSE, etINT,{&nmultisim},
3610 "HIDDENDo multiple simulations in parallel" },
3611 { "-replex", FALSE, etINT, {&repl_ex_nst},
3612 "HIDDENAttempt replica exchange periodically with this period (steps)" },
3613 { "-reseed", FALSE, etINT, {&repl_ex_seed},
3614 "HIDDENSeed for replica exchange, -1 is generate a seed" },
3615 { "-rerunvsite", FALSE, etBOOL, {&bRerunVSite},
3616 "HIDDENRecalculate virtual site coordinates with [TT]-rerun[tt]" },
3617 { "-ionize", FALSE, etBOOL,{&bIonize},
3618 "HIDDENDo a simulation including the effect of an X-Ray bombardment on your system" },
3619 { "-confout", TRUE, etBOOL, {&bConfout},
3620 "HIDDENWrite the last configuration with [TT]-c[tt] and force checkpointing at the last step" },
3621 { "-stepout", FALSE, etINT, {&nstepout},
3622 "HIDDENFrequency of writing the remaining runtime" },
3623 { "-resetstep", FALSE, etINT, {&resetstep},
3624 "HIDDENReset cycle counters after these many time steps" },
3625 { "-resethway", FALSE, etBOOL, {&bResetCountersHalfWay},
3626 "HIDDENReset the cycle counters after half the number of steps or halfway [TT]-maxh[tt]" },
3627 { "-v", FALSE, etBOOL,{&bVerbose},
3628 "Be loud and noisy" },
3629 { "-maxh", FALSE, etREAL, {&max_hours},
3630 "HIDDENTerminate after 0.99 times this time (hours)" },
3631 { "-cpt", FALSE, etREAL, {&cpt_period},
3632 "HIDDENCheckpoint interval (minutes)" },
3633 { "-append", FALSE, etBOOL, {&bAppendFiles},
3634 "HIDDENAppend to previous output files when continuing from checkpoint" },
3635 { "-addpart", FALSE, etBOOL, {&bAddPart},
3636 "HIDDENAdd the simulation part number to all output files when continuing from checkpoint" },
3637 };
3638 gmx_edsam_t ed;
3639 unsigned long Flags, PCA_Flags;
3640 ivec ddxyz;
3641 int dd_node_order;
3642 gmx_bool HaveCheckpoint;
3643 FILE *fplog,*fptest;
3644 int sim_part,sim_part_fn;
3645 const char *part_suffix=".part";
3646 char suffix[STRLEN];
3647 int rc;
3648 char **multidir=NULL;
3649
3650 cr = init_par(&argc,&argv);
3651
3652 PCA_Flags = (PCA_KEEP_ARGS | PCA_NOEXIT_ON_ARGS | PCA_CAN_SET_DEFFNM
3653 | (MASTER(cr) ? 0 : PCA_QUIET));
3654
3655
3656 /* Comment this in to do fexist calls only on master
3657 * works not with rerun or tables at the moment
3658 * also comment out the version of init_forcerec in md.c
3659 * with NULL instead of opt2fn
3660 */
3661 /*
3662 if (!MASTER(cr))
3663 {
3664 PCA_Flags |= PCA_NOT_READ_NODE;
3665 }
3666 */
3667
3668 parse_common_args(&argc,argv,PCA_Flags, NFILE,fnm,asize(pa),pa,
3669 asize(desc),desc,0,NULL, &oenv);
3670
3671 /* we set these early because they might be used in init_multisystem()
3672 Note that there is the potential for npme>nnodes until the number of
3673 threads is set later on, if there's thread parallelization. That shouldn't
3674 lead to problems. */
3675 dd_node_order = nenum(ddno_opt);
3676 cr->npmenodes = npme;
3677
3678 #ifdef GMX_THREADS
3679 /* now determine the number of threads automatically. The threads are
3680 only started at mdrunner_threads, though. */
3681 if (nthreads<1)
3682 {
3683 nthreads=tMPI_Thread_get_hw_number();
3684 }
3685 #else
3686 nthreads=1;
3687 #endif
3688
3689 /* now check the -multi and -multidir option */
3690 if (opt2bSet("-multidir", NFILE, fnm))
3691 {
3692 int i;
3693 if (nmultisim > 0)
3694 {
3695 gmx_fatal(FARGS, "mdrun -multi and -multidir options are mutually exclusive.");
3696 }
3697 nmultisim = opt2fns(&multidir, "-multidir", NFILE, fnm);
3698 }
3699
3700
3701 if (repl_ex_nst != 0 && nmultisim < 2)
3702 gmx_fatal(FARGS,"Need at least two replicas for replica exchange (option -multi)");
3703
3704 if (nmultisim > 1) {
3705 #ifndef GMX_THREADS
3706 gmx_bool bParFn = (multidir == NULL);
3707 init_multisystem(cr,nmultisim,multidir,NFILE,fnm,TRUE);
3708 #else
3709 gmx_fatal(FARGS,"mdrun -multi is not supported with the thread library.Please compile GROMACS with MPI support");
3710 #endif
3711 }
3712
3713 /* Check if there is ANY checkpoint file available */
3714 sim_part = 1;
3715 sim_part_fn = sim_part;
3716 if (opt2bSet("-cpi",NFILE,fnm))
3717 {
3718 bAppendFiles =
3719 read_checkpoint_simulation_part(opt2fn_master("-cpi", NFILE,fnm,cr),
3720 &sim_part_fn,NULL,cr,
3721 bAppendFiles,NFILE,fnm,
3722 part_suffix,&bAddPart);
3723 if (sim_part_fn==0 && MASTER(cr))
3724 {
3725 fprintf(stdout,"No previous checkpoint file present, assuming this is a new run.\n");
3726 }
3727 else
3728 {
3729 sim_part = sim_part_fn + 1;
3730 }
3731 }
3732 else
3733 {
3734 bAppendFiles = FALSE;
3735 }
3736
3737 if (!bAppendFiles)
3738 {
3739 sim_part_fn = sim_part;
3740 }
3741
3742 if (bAddPart && sim_part_fn > 1)
3743 {
3744 /* This is a continuation run, rename trajectory output files
3745 (except checkpoint files) */
3746 /* create new part name first (zero-filled) */
3747 sprintf(suffix,"%s%04d",part_suffix,sim_part_fn);
3748
3749 add_suffix_to_output_names(fnm,NFILE,suffix);
3750 fprintf(stdout,"Checkpoint file is from part %d, new output files will be suffixed '%s'.\n",sim_part-1,suffix);
3751 }
3752
3753 Flags = opt2bSet("-rerun",NFILE,fnm) ? MD_RERUN : 0;
3754 Flags = Flags | (bSepPot ? MD_SEPPOT : 0);
3755 Flags = Flags | (bIonize ? MD_IONIZE : 0);
3756 Flags = Flags | (bPartDec ? MD_PARTDEC : 0);
3757 Flags = Flags | (bDDBondCheck ? MD_DDBONDCHECK : 0);
3758 Flags = Flags | (bDDBondComm ? MD_DDBONDCOMM : 0);
3759 Flags = Flags | (bConfout ? MD_CONFOUT : 0);
3760 Flags = Flags | (bRerunVSite ? MD_RERUN_VSITE : 0);
3761 Flags = Flags | (bReproducible ? MD_REPRODUCIBLE : 0);
3762 Flags = Flags | (bAppendFiles ? MD_APPENDFILES : 0);
3763 Flags = Flags | (opt2parg_bSet("-append", asize(pa),pa) ? MD_APPENDFILESSET : 0);
3764 Flags = Flags | (sim_part>1 ? MD_STARTFROMCPT : 0);
3765 Flags = Flags | (bResetCountersHalfWay ? MD_RESETCOUNTERSHALFWAY : 0);
3766
3767
3768 /* We postpone opening the log file if we are appending, so we can
3769 first truncate the old log file and append to the correct position
3770 there instead. */
3771 if ((MASTER(cr) || bSepPot) && !bAppendFiles)
3772 {
3773 gmx_log_open(ftp2fn(efLOG,NFILE,fnm),cr,!bSepPot,Flags,&fplog);
3774 CopyRight(fplog,argv[0]);
3775 please_cite(fplog,"Hess2008b");
3776 please_cite(fplog,"Spoel2005a");
3777 please_cite(fplog,"Lindahl2001a");
3778 please_cite(fplog,"Berendsen95a");
3779 }
3780 else
3781 {
3782 fplog = NULL;
3783 }
3784
3785 ddxyz[XX] = (int)(realddxyz[XX] + 0.5);
3786 ddxyz[YY] = (int)(realddxyz[YY] + 0.5);
3787 ddxyz[ZZ] = (int)(realddxyz[ZZ] + 0.5);
3788
3789 /* even if nthreads = 1, we still call this one */
3790
3791 rc = mdrunner_membed(fplog, cr, NFILE, fnm, oenv, bVerbose, bCompact,
3792 nstglobalcomm,
3793 ddxyz, dd_node_order, rdd, rconstr, dddlb_opt[0], dlb_scale,
3794 ddcsx, ddcsy, ddcsz, nstepout, resetstep, nmultisim, repl_ex_nst,
3795 repl_ex_seed, pforce, cpt_period, max_hours, deviceOptions, Flags,
3796 xy_fac,xy_max,z_fac,z_max,
3797 it_xy,it_z,probe_rad,low_up_rm,
3798 pieces,bALLOW_ASYMMETRY,maxwarn);
3799
3800 if (gmx_parallel_env_initialized())
3801 gmx_finalize();
3802
3803 if (MULTIMASTER(cr)) {
3804 thanx(stderr);
3805 }
3806
3807 /* Log file has to be closed in mdrunner if we are appending to it
3808 (fplog not set here) */
3809 fprintf(stderr,"Please cite:\nWolf et al, J Comp Chem 31 (2010) 2169-2174.\n");
3810
3811 if (MASTER(cr) && !bAppendFiles)
3812 {
3813 gmx_log_close(fplog);
3814 }
3815
3816 return rc;
3817 }
The ultimate molecular dynamics simulation package