3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
41 #include "gmx_fatal.h"
56 #include "mtop_util.h"
59 #define RANGECHK(i,n) if ((i)>=(n)) gmx_fatal(FARGS,"Your index file contains atomnumbers (e.g. %d)\nthat are larger than the number of atoms in the tpr file (%d)",(i),(n))
61 static gmx_bool *bKeepIt(int gnx,int natoms,atom_id index[])
67 for(i=0; (i<gnx); i++) {
68 RANGECHK(index[i],natoms);
75 static atom_id *invind(int gnx,int natoms,atom_id index[])
81 for(i=0; (i<gnx); i++) {
82 RANGECHK(index[i],natoms);
89 static void reduce_block(gmx_bool bKeep[],t_block *block,
95 snew(index,block->nr);
98 for(i=0; (i<block->nr); i++) {
99 for(j=block->index[i]; (j<block->index[i+1]); j++) {
104 if (newj > index[newi]) {
110 fprintf(stderr,"Reduced block %8s from %6d to %6d index-, %6d to %6d a-entries\n",
111 name,block->nr,newi,block->index[block->nr],newj);
112 block->index = index;
116 static void reduce_blocka(atom_id invindex[],gmx_bool bKeep[],t_blocka *block,
122 snew(index,block->nr);
126 for(i=0; (i<block->nr); i++) {
127 for(j=block->index[i]; (j<block->index[i+1]); j++) {
130 a[newj] = invindex[k];
134 if (newj > index[newi]) {
140 fprintf(stderr,"Reduced block %8s from %6d to %6d index-, %6d to %6d a-entries\n",
141 name,block->nr,newi,block->nra,newj);
142 block->index = index;
148 static void reduce_rvec(int gnx,atom_id index[],rvec vv[])
154 for(i=0; (i<gnx); i++)
155 copy_rvec(vv[index[i]],ptr[i]);
156 for(i=0; (i<gnx); i++)
157 copy_rvec(ptr[i],vv[i]);
161 static void reduce_atom(int gnx,atom_id index[],t_atom atom[],char ***atomname,
162 int *nres,t_resinfo *resinfo)
171 snew(rinfo,atom[index[gnx-1]].resind+1);
172 for(i=0; (i<gnx); i++) {
173 ptr[i] = atom[index[i]];
174 aname[i] = atomname[index[i]];
177 for(i=0; (i<gnx); i++) {
179 atomname[i] = aname[i];
180 if ((i==0) || (atom[i].resind != atom[i-1].resind)) {
182 rinfo[nr] = resinfo[atom[i].resind];
187 for(i=0; (i<nr); i++) {
188 resinfo[i] = rinfo[i];
197 static void reduce_ilist(atom_id invindex[],gmx_bool bKeep[],
198 t_ilist *il,int nratoms,const char *name)
207 for(i=0; (i<il->nr); i+=nratoms+1) {
209 for(j=1; (j<=nratoms); j++) {
210 bB = bB && bKeep[il->iatoms[i+j]];
213 ia[newnr++] = il->iatoms[i];
214 for(j=1; (j<=nratoms); j++)
215 ia[newnr++] = invindex[il->iatoms[i+j]];
218 fprintf(stderr,"Reduced ilist %8s from %6d to %6d entries\n",
219 name,il->nr/(nratoms+1),
223 for(i=0; (i<newnr); i++)
224 il->iatoms[i] = ia[i];
230 static void reduce_topology_x(int gnx,atom_id index[],
231 gmx_mtop_t *mtop,rvec x[],rvec v[])
238 top = gmx_mtop_t_to_t_topology(mtop);
239 bKeep = bKeepIt(gnx,top.atoms.nr,index);
240 invindex = invind(gnx,top.atoms.nr,index);
242 reduce_block(bKeep,&(top.cgs),"cgs");
243 reduce_block(bKeep,&(top.mols),"mols");
244 reduce_blocka(invindex,bKeep,&(top.excls),"excls");
245 reduce_rvec(gnx,index,x);
246 reduce_rvec(gnx,index,v);
247 reduce_atom(gnx,index,top.atoms.atom,top.atoms.atomname,
248 &(top.atoms.nres),top.atoms.resinfo);
250 for(i=0; (i<F_NRE); i++) {
251 reduce_ilist(invindex,bKeep,&(top.idef.il[i]),
252 interaction_function[i].nratoms,
253 interaction_function[i].name);
259 snew(mtop->moltype,mtop->nmoltype);
260 mtop->moltype[0].name = mtop->name;
261 mtop->moltype[0].atoms = top.atoms;
262 for(i=0; i<F_NRE; i++) {
263 mtop->moltype[0].ilist[i] = top.idef.il[i];
265 mtop->moltype[0].atoms = top.atoms;
266 mtop->moltype[0].cgs = top.cgs;
267 mtop->moltype[0].excls = top.excls;
270 snew(mtop->molblock,mtop->nmolblock);
271 mtop->molblock[0].type = 0;
272 mtop->molblock[0].nmol = 1;
273 mtop->molblock[0].natoms_mol = top.atoms.nr;
274 mtop->molblock[0].nposres_xA = 0;
275 mtop->molblock[0].nposres_xB = 0;
277 mtop->natoms = top.atoms.nr;
280 static void zeroq(int n,atom_id index[],gmx_mtop_t *mtop)
284 for(mt=0; mt<mtop->nmoltype; mt++) {
285 for(i=0; (i<mtop->moltype[mt].atoms.nr); i++) {
286 mtop->moltype[mt].atoms.atom[index[i]].q = 0;
287 mtop->moltype[mt].atoms.atom[index[i]].qB = 0;
292 int main (int argc, char *argv[])
294 const char *desc[] = {
295 "tpbconv can edit run input files in four ways.[PAR]",
296 "[BB]1st.[bb] by modifying the number of steps in a run input file",
297 "with options [TT]-extend[tt], [TT]-until[tt] or [TT]-nsteps[tt]",
298 "(nsteps=-1 means unlimited number of steps)[PAR]",
299 "[BB]2nd.[bb] (OBSOLETE) by creating a run input file",
300 "for a continuation run when your simulation has crashed due to e.g.",
301 "a full disk, or by making a continuation run input file.",
302 "This option is obsolete, since mdrun now writes and reads",
304 "Note that a frame with coordinates and velocities is needed.",
305 "When pressure and/or Nose-Hoover temperature coupling is used",
306 "an energy file can be supplied to get an exact continuation",
307 "of the original run.[PAR]",
308 "[BB]3rd.[bb] by creating a tpx file for a subset of your original",
309 "tpx file, which is useful when you want to remove the solvent from",
310 "your tpx file, or when you want to make e.g. a pure Ca tpx file.",
311 "Note that you may need to use -nsteps -1 (or similar) to get",
313 "[BB]WARNING: this tpx file is not fully functional[bb].[PAR]",
314 "[BB]4th.[bb] by setting the charges of a specified group",
315 "to zero. This is useful when doing free energy estimates",
316 "using the LIE (Linear Interaction Energy) method."
319 const char *top_fn,*frame_fn;
321 ener_file_t fp_ener=NULL;
324 gmx_large_int_t nsteps_req,run_step,frame;
325 double run_t,state_t;
326 gmx_bool bOK,bNsteps,bExtend,bUntil,bTime,bTraj;
327 gmx_bool bFrame,bUse,bSel,bNeedEner,bReadEner,bScanEner;
330 t_inputrec *ir,*irnew=NULL;
333 rvec *newx=NULL,*newv=NULL,*tmpx,*tmpv;
339 gmx_enxnm_t *enm=NULL;
340 t_enxframe *fr_ener=NULL;
341 char buf[200],buf2[200];
344 { efTPX, NULL, NULL, ffREAD },
345 { efTRN, "-f", NULL, ffOPTRD },
346 { efEDR, "-e", NULL, ffOPTRD },
347 { efNDX, NULL, NULL, ffOPTRD },
348 { efTPX, "-o", "tpxout",ffWRITE }
350 #define NFILE asize(fnm)
352 /* Command line options */
353 static int nsteps_req_int = 0;
354 static real start_t = -1.0, extend_t = 0.0, until_t = 0.0;
355 static gmx_bool bContinuation = TRUE,bZeroQ = FALSE,bVel=TRUE;
356 static t_pargs pa[] = {
357 { "-extend", FALSE, etREAL, {&extend_t},
358 "Extend runtime by this amount (ps)" },
359 { "-until", FALSE, etREAL, {&until_t},
360 "Extend runtime until this ending time (ps)" },
361 { "-nsteps", FALSE, etINT, {&nsteps_req_int},
362 "Change the number of steps" },
363 { "-time", FALSE, etREAL, {&start_t},
364 "Continue from frame at this time (ps) instead of the last frame" },
365 { "-zeroq", FALSE, etBOOL, {&bZeroQ},
366 "Set the charges of a group (from the index) to zero" },
367 { "-vel", FALSE, etBOOL, {&bVel},
368 "Require velocities from trajectory" },
369 { "-cont", FALSE, etBOOL, {&bContinuation},
370 "For exact continuation, the constraints should not be applied before the first step" }
374 CopyRight(stdout,argv[0]);
376 /* Parse the command line */
377 parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
378 asize(desc),desc,0,NULL,&oenv);
380 /* Convert int to gmx_large_int_t */
381 nsteps_req = nsteps_req_int;
382 bNsteps = opt2parg_bSet("-nsteps",asize(pa),pa);
383 bExtend = opt2parg_bSet("-extend",asize(pa),pa);
384 bUntil = opt2parg_bSet("-until",asize(pa),pa);
385 bTime = opt2parg_bSet("-time",asize(pa),pa);
386 bTraj = (opt2bSet("-f",NFILE,fnm) || bTime);
388 top_fn = ftp2fn(efTPX,NFILE,fnm);
389 fprintf(stderr,"Reading toplogy and stuff from %s\n",top_fn);
392 read_tpx_state(top_fn,ir,&state,NULL,&mtop);
393 run_step = ir->init_step;
394 run_t = ir->init_step*ir->delta_t + ir->init_t;
396 if (!EI_STATE_VELOCITY(ir->eI)) {
402 "NOTE: Reading the state from trajectory is an obsolete feaure of tpbconv.\n"
403 " Continuation should be done by loading a checkpoint file with mdrun -cpi\n"
404 " This guarantees that all state variables are transferred.\n"
405 " tpbconv is now only useful for increasing nsteps,\n"
406 " but even that can often be avoided by using mdrun -maxh\n"
409 if (ir->bContinuation != bContinuation)
410 fprintf(stderr,"Modifying ir->bContinuation to %s\n",
411 bool_names[bContinuation]);
412 ir->bContinuation = bContinuation;
415 bNeedEner = (ir->epc == epcPARRINELLORAHMAN || ir->etc == etcNOSEHOOVER);
416 bReadEner = (bNeedEner && ftp2bSet(efEDR,NFILE,fnm));
417 bScanEner = (bReadEner && !bTime);
419 if (ir->epc != epcNO || EI_SD(ir->eI) || ir->eI == eiBD) {
420 fprintf(stderr,"NOTE: The simulation uses pressure coupling and/or stochastic dynamics.\n"
421 "tpbconv can not provide binary identical continuation.\n"
422 "If you want that, supply a checkpoint file to mdrun\n\n");
425 if (EI_SD(ir->eI) || ir->eI == eiBD) {
426 fprintf(stderr,"\nChanging ld-seed from %d ",ir->ld_seed);
427 ir->ld_seed = make_seed();
428 fprintf(stderr,"to %d\n\n",ir->ld_seed);
431 frame_fn = ftp2fn(efTRN,NFILE,fnm);
433 "\nREADING COORDS, VELS AND BOX FROM TRAJECTORY %s...\n\n",
436 fp = open_trn(frame_fn,"r");
438 fp_ener = open_enx(ftp2fn(efEDR,NFILE,fnm),"r");
439 do_enxnms(fp_ener,&nre,&enm);
444 /* Now scan until the last set of x and v (step == 0)
445 * or the ones at step step.
450 bFrame = fread_trnheader(fp,&head,&bOK);
451 if (bOK && frame == 0) {
452 if (mtop.natoms != head.natoms)
453 gmx_fatal(FARGS,"Number of atoms in Topology (%d) "
454 "is not the same as in Trajectory (%d)\n",
455 mtop.natoms,head.natoms);
456 snew(newx,head.natoms);
457 snew(newv,head.natoms);
459 bFrame = bFrame && bOK;
462 bOK = fread_htrn(fp,&head,newbox,newx,newv,NULL);
464 bFrame = bFrame && bOK;
467 (head.x_size) && (head.v_size || !bVel)) {
470 /* Read until the energy time is >= the trajectory time */
471 while (fr_ener->t < head.t && do_enx(fp_ener,fr_ener));
472 bUse = (fr_ener->t == head.t);
482 run_step = head.step;
483 state.lambda = head.lambda;
484 copy_mat(newbox,state.box);
487 if (bFrame || !bOK) {
488 sprintf(buf,"\r%s %s frame %s%s: step %s%s time %s",
489 "%s","%s","%6",gmx_large_int_fmt,"%6",gmx_large_int_fmt," %8.3f");
491 bUse ? "Read " : "Skipped",ftp2ext(fn2ftp(frame_fn)),
492 frame,head.step,head.t);
494 if (bTime && (head.t >= start_t))
500 free_enxframe(fr_ener);
501 free_enxnms(nre,enm);
504 fprintf(stderr,"\n");
507 fprintf(stderr,"%s frame %s (step %s, time %g) is incomplete\n",
508 ftp2ext(fn2ftp(frame_fn)),gmx_step_str(frame-1,buf2),
509 gmx_step_str(head.step,buf),head.t);
510 fprintf(stderr,"\nUsing frame of step %s time %g\n",
511 gmx_step_str(run_step,buf),run_t);
515 get_enx_state(ftp2fn(efEDR,NFILE,fnm),run_t,&mtop.groups,ir,&state);
517 fprintf(stderr,"\nWARNING: The simulation uses %s temperature and/or %s pressure coupling,\n"
518 " the continuation will only be exact when an energy file is supplied\n\n",
519 ETCOUPLTYPE(etcNOSEHOOVER),
520 EPCOUPLTYPE(epcPARRINELLORAHMAN));
526 fprintf(stderr,"Setting nsteps to %s\n",gmx_step_str(nsteps_req,buf));
527 ir->nsteps = nsteps_req;
529 /* Determine total number of steps remaining */
531 ir->nsteps = ir->nsteps - (run_step - ir->init_step) + (gmx_large_int_t)(extend_t/ir->delta_t + 0.5);
532 printf("Extending remaining runtime of by %g ps (now %s steps)\n",
533 extend_t,gmx_step_str(ir->nsteps,buf));
536 printf("nsteps = %s, run_step = %s, current_t = %g, until = %g\n",
537 gmx_step_str(ir->nsteps,buf),
538 gmx_step_str(run_step,buf2),
540 ir->nsteps = (gmx_large_int_t)((until_t - run_t)/ir->delta_t + 0.5);
541 printf("Extending remaining runtime until %g ps (now %s steps)\n",
542 until_t,gmx_step_str(ir->nsteps,buf));
545 ir->nsteps -= run_step - ir->init_step;
547 printf("%s steps (%g ps) remaining from first run.\n",
548 gmx_step_str(ir->nsteps,buf),ir->nsteps*ir->delta_t);
553 if (bNsteps || bZeroQ || (ir->nsteps > 0)) {
554 ir->init_step = run_step;
556 if (ftp2bSet(efNDX,NFILE,fnm) ||
557 !(bNsteps || bExtend || bUntil || bTraj)) {
558 atoms = gmx_mtop_global_atoms(&mtop);
559 get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),1,
560 &gnx,&index,&grpname);
562 bSel = (gnx != state.natoms);
563 for (i=0; ((i<gnx) && (!bSel)); i++)
564 bSel = (i!=index[i]);
569 fprintf(stderr,"Will write subset %s of original tpx containing %d "
570 "atoms\n",grpname,gnx);
571 reduce_topology_x(gnx,index,&mtop,state.x,state.v);
575 zeroq(gnx,index,&mtop);
576 fprintf(stderr,"Zero-ing charges for group %s\n",grpname);
579 fprintf(stderr,"Will write full tpx file (no selection)\n");
582 state_t = ir->init_t + ir->init_step*ir->delta_t;
583 sprintf(buf, "Writing statusfile with starting step %s%s and length %s%s steps...\n","%10",gmx_large_int_fmt,"%10",gmx_large_int_fmt);
584 fprintf(stderr,buf,ir->init_step,ir->nsteps);
585 fprintf(stderr," time %10.3f and length %10.3f ps\n",
586 state_t,ir->nsteps*ir->delta_t);
587 write_tpx_state(opt2fn("-o",NFILE,fnm),ir,&state,&mtop);
590 printf("You've simulated long enough. Not writing tpr file\n");