2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
41 #include "gromacs/legacyheaders/types/inputrec.h"
42 #include "gromacs/legacyheaders/types/simple.h"
43 #include "gromacs/legacyheaders/types/state.h"
44 #include "gromacs/topology/index.h"
47 #include "gromacs/gmxpreprocess/readir.h"
48 #include "gromacs/topology/mtop_util.h"
49 #include "checkpoint.h"
50 #include "gromacs/fileio/tpxio.h"
51 #include "gromacs/fileio/trnio.h"
52 #include "gromacs/fileio/enxio.h"
53 #include "gromacs/utility/futil.h"
55 #include "gromacs/commandline/pargs.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/random/random.h"
58 #include "gromacs/topology/topology.h"
59 #include "gromacs/utility/fatalerror.h"
60 #include "gromacs/utility/smalloc.h"
62 #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))
64 static gmx_bool
*bKeepIt(int gnx
, int natoms
, atom_id index
[])
70 for (i
= 0; (i
< gnx
); i
++)
72 RANGECHK(index
[i
], natoms
);
79 static atom_id
*invind(int gnx
, int natoms
, atom_id index
[])
85 for (i
= 0; (i
< gnx
); i
++)
87 RANGECHK(index
[i
], natoms
);
94 static void reduce_block(gmx_bool bKeep
[], t_block
*block
,
100 snew(index
, block
->nr
);
103 for (i
= 0; (i
< block
->nr
); i
++)
105 for (j
= block
->index
[i
]; (j
< block
->index
[i
+1]); j
++)
112 if (newj
> index
[newi
])
119 fprintf(stderr
, "Reduced block %8s from %6d to %6d index-, %6d to %6d a-entries\n",
120 name
, block
->nr
, newi
, block
->index
[block
->nr
], newj
);
121 block
->index
= index
;
125 static void reduce_blocka(atom_id invindex
[], gmx_bool bKeep
[], t_blocka
*block
,
129 int i
, j
, k
, newi
, newj
;
131 snew(index
, block
->nr
);
135 for (i
= 0; (i
< block
->nr
); i
++)
137 for (j
= block
->index
[i
]; (j
< block
->index
[i
+1]); j
++)
142 a
[newj
] = invindex
[k
];
146 if (newj
> index
[newi
])
153 fprintf(stderr
, "Reduced block %8s from %6d to %6d index-, %6d to %6d a-entries\n",
154 name
, block
->nr
, newi
, block
->nra
, newj
);
155 block
->index
= index
;
161 static void reduce_rvec(int gnx
, atom_id index
[], rvec vv
[])
167 for (i
= 0; (i
< gnx
); i
++)
169 copy_rvec(vv
[index
[i
]], ptr
[i
]);
171 for (i
= 0; (i
< gnx
); i
++)
173 copy_rvec(ptr
[i
], vv
[i
]);
178 static void reduce_atom(int gnx
, atom_id index
[], t_atom atom
[], char ***atomname
,
179 int *nres
, t_resinfo
*resinfo
)
188 snew(rinfo
, atom
[index
[gnx
-1]].resind
+1);
189 for (i
= 0; (i
< gnx
); i
++)
191 ptr
[i
] = atom
[index
[i
]];
192 aname
[i
] = atomname
[index
[i
]];
195 for (i
= 0; (i
< gnx
); i
++)
198 atomname
[i
] = aname
[i
];
199 if ((i
== 0) || (atom
[i
].resind
!= atom
[i
-1].resind
))
202 rinfo
[nr
] = resinfo
[atom
[i
].resind
];
207 for (i
= 0; (i
< nr
); i
++)
209 resinfo
[i
] = rinfo
[i
];
218 static void reduce_ilist(atom_id invindex
[], gmx_bool bKeep
[],
219 t_ilist
*il
, int nratoms
, const char *name
)
229 for (i
= 0; (i
< il
->nr
); i
+= nratoms
+1)
232 for (j
= 1; (j
<= nratoms
); j
++)
234 bB
= bB
&& bKeep
[il
->iatoms
[i
+j
]];
238 ia
[newnr
++] = il
->iatoms
[i
];
239 for (j
= 1; (j
<= nratoms
); j
++)
241 ia
[newnr
++] = invindex
[il
->iatoms
[i
+j
]];
245 fprintf(stderr
, "Reduced ilist %8s from %6d to %6d entries\n",
246 name
, il
->nr
/(nratoms
+1),
250 for (i
= 0; (i
< newnr
); i
++)
252 il
->iatoms
[i
] = ia
[i
];
259 static void reduce_topology_x(int gnx
, atom_id index
[],
260 gmx_mtop_t
*mtop
, rvec x
[], rvec v
[])
267 top
= gmx_mtop_t_to_t_topology(mtop
);
268 bKeep
= bKeepIt(gnx
, top
.atoms
.nr
, index
);
269 invindex
= invind(gnx
, top
.atoms
.nr
, index
);
271 reduce_block(bKeep
, &(top
.cgs
), "cgs");
272 reduce_block(bKeep
, &(top
.mols
), "mols");
273 reduce_blocka(invindex
, bKeep
, &(top
.excls
), "excls");
274 reduce_rvec(gnx
, index
, x
);
275 reduce_rvec(gnx
, index
, v
);
276 reduce_atom(gnx
, index
, top
.atoms
.atom
, top
.atoms
.atomname
,
277 &(top
.atoms
.nres
), top
.atoms
.resinfo
);
279 for (i
= 0; (i
< F_NRE
); i
++)
281 reduce_ilist(invindex
, bKeep
, &(top
.idef
.il
[i
]),
282 interaction_function
[i
].nratoms
,
283 interaction_function
[i
].name
);
289 snew(mtop
->moltype
, mtop
->nmoltype
);
290 mtop
->moltype
[0].name
= mtop
->name
;
291 mtop
->moltype
[0].atoms
= top
.atoms
;
292 for (i
= 0; i
< F_NRE
; i
++)
294 mtop
->moltype
[0].ilist
[i
] = top
.idef
.il
[i
];
296 mtop
->moltype
[0].atoms
= top
.atoms
;
297 mtop
->moltype
[0].cgs
= top
.cgs
;
298 mtop
->moltype
[0].excls
= top
.excls
;
301 snew(mtop
->molblock
, mtop
->nmolblock
);
302 mtop
->molblock
[0].type
= 0;
303 mtop
->molblock
[0].nmol
= 1;
304 mtop
->molblock
[0].natoms_mol
= top
.atoms
.nr
;
305 mtop
->molblock
[0].nposres_xA
= 0;
306 mtop
->molblock
[0].nposres_xB
= 0;
308 mtop
->natoms
= top
.atoms
.nr
;
311 static void zeroq(atom_id index
[], gmx_mtop_t
*mtop
)
315 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
317 for (i
= 0; (i
< mtop
->moltype
[mt
].atoms
.nr
); i
++)
319 mtop
->moltype
[mt
].atoms
.atom
[index
[i
]].q
= 0;
320 mtop
->moltype
[mt
].atoms
.atom
[index
[i
]].qB
= 0;
325 int gmx_convert_tpr(int argc
, char *argv
[])
327 const char *desc
[] = {
328 "[THISMODULE] can edit run input files in four ways.[PAR]",
329 "[BB]1.[bb] by modifying the number of steps in a run input file",
330 "with options [TT]-extend[tt], [TT]-until[tt] or [TT]-nsteps[tt]",
331 "(nsteps=-1 means unlimited number of steps)[PAR]",
332 "[BB]2.[bb] (OBSOLETE) by creating a run input file",
333 "for a continuation run when your simulation has crashed due to e.g.",
334 "a full disk, or by making a continuation run input file.",
335 "This option is obsolete, since mdrun now writes and reads",
337 "[BB]Note[bb] that a frame with coordinates and velocities is needed.",
338 "When pressure and/or Nose-Hoover temperature coupling is used",
339 "an energy file can be supplied to get an exact continuation",
340 "of the original run.[PAR]",
341 "[BB]3.[bb] by creating a [TT].tpx[tt] file for a subset of your original",
342 "tpx file, which is useful when you want to remove the solvent from",
343 "your [TT].tpx[tt] file, or when you want to make e.g. a pure C[GRK]alpha[grk] [TT].tpx[tt] file.",
344 "Note that you may need to use [TT]-nsteps -1[tt] (or similar) to get",
346 "[BB]WARNING: this [TT].tpx[tt] file is not fully functional[bb].[PAR]",
347 "[BB]4.[bb] by setting the charges of a specified group",
348 "to zero. This is useful when doing free energy estimates",
349 "using the LIE (Linear Interaction Energy) method."
352 const char *top_fn
, *frame_fn
;
354 ener_file_t fp_ener
= NULL
;
357 gmx_int64_t nsteps_req
, run_step
, frame
;
358 double run_t
, state_t
;
359 gmx_bool bOK
, bNsteps
, bExtend
, bUntil
, bTime
, bTraj
;
360 gmx_bool bFrame
, bUse
, bSel
, bNeedEner
, bReadEner
, bScanEner
, bFepState
;
363 t_inputrec
*ir
, *irnew
= NULL
;
366 rvec
*newx
= NULL
, *newv
= NULL
, *tmpx
, *tmpv
;
370 atom_id
*index
= NULL
;
372 gmx_enxnm_t
*enm
= NULL
;
373 t_enxframe
*fr_ener
= NULL
;
374 char buf
[200], buf2
[200];
377 { efTPX
, NULL
, NULL
, ffREAD
},
378 { efTRN
, "-f", NULL
, ffOPTRD
},
379 { efEDR
, "-e", NULL
, ffOPTRD
},
380 { efNDX
, NULL
, NULL
, ffOPTRD
},
381 { efTPX
, "-o", "tpxout", ffWRITE
}
383 #define NFILE asize(fnm)
385 /* Command line options */
386 static int nsteps_req_int
= 0;
387 static real start_t
= -1.0, extend_t
= 0.0, until_t
= 0.0;
388 static int init_fep_state
= 0;
389 static gmx_bool bContinuation
= TRUE
, bZeroQ
= FALSE
, bVel
= TRUE
;
390 static t_pargs pa
[] = {
391 { "-extend", FALSE
, etREAL
, {&extend_t
},
392 "Extend runtime by this amount (ps)" },
393 { "-until", FALSE
, etREAL
, {&until_t
},
394 "Extend runtime until this ending time (ps)" },
395 { "-nsteps", FALSE
, etINT
, {&nsteps_req_int
},
396 "Change the number of steps" },
397 { "-time", FALSE
, etREAL
, {&start_t
},
398 "Continue from frame at this time (ps) instead of the last frame" },
399 { "-zeroq", FALSE
, etBOOL
, {&bZeroQ
},
400 "Set the charges of a group (from the index) to zero" },
401 { "-vel", FALSE
, etBOOL
, {&bVel
},
402 "Require velocities from trajectory" },
403 { "-cont", FALSE
, etBOOL
, {&bContinuation
},
404 "For exact continuation, the constraints should not be applied before the first step" },
405 { "-init_fep_state", FALSE
, etINT
, {&init_fep_state
},
406 "fep state to initialize from" },
410 /* Parse the command line */
411 if (!parse_common_args(&argc
, argv
, 0, NFILE
, fnm
, asize(pa
), pa
,
412 asize(desc
), desc
, 0, NULL
, &oenv
))
417 /* Convert int to gmx_int64_t */
418 nsteps_req
= nsteps_req_int
;
419 bNsteps
= opt2parg_bSet("-nsteps", asize(pa
), pa
);
420 bExtend
= opt2parg_bSet("-extend", asize(pa
), pa
);
421 bUntil
= opt2parg_bSet("-until", asize(pa
), pa
);
422 bFepState
= opt2parg_bSet("-init_fep_state", asize(pa
), pa
);
423 bTime
= opt2parg_bSet("-time", asize(pa
), pa
);
424 bTraj
= (opt2bSet("-f", NFILE
, fnm
) || bTime
);
426 top_fn
= ftp2fn(efTPX
, NFILE
, fnm
);
427 fprintf(stderr
, "Reading toplogy and stuff from %s\n", top_fn
);
430 read_tpx_state(top_fn
, ir
, &state
, NULL
, &mtop
);
431 run_step
= ir
->init_step
;
432 run_t
= ir
->init_step
*ir
->delta_t
+ ir
->init_t
;
434 if (!EI_STATE_VELOCITY(ir
->eI
))
442 "NOTE: Reading the state from trajectory is an obsolete feature of gmx convert-tpr.\n"
443 " Continuation should be done by loading a checkpoint file with mdrun -cpi\n"
444 " This guarantees that all state variables are transferred.\n"
445 " gmx convert-tpr is now only useful for increasing nsteps,\n"
446 " but even that can often be avoided by using mdrun -maxh\n"
449 if (ir
->bContinuation
!= bContinuation
)
451 fprintf(stderr
, "Modifying ir->bContinuation to %s\n",
452 bool_names
[bContinuation
]);
454 ir
->bContinuation
= bContinuation
;
457 bNeedEner
= (ir
->epc
== epcPARRINELLORAHMAN
|| ir
->etc
== etcNOSEHOOVER
);
458 bReadEner
= (bNeedEner
&& ftp2bSet(efEDR
, NFILE
, fnm
));
459 bScanEner
= (bReadEner
&& !bTime
);
461 if (ir
->epc
!= epcNO
|| EI_SD(ir
->eI
) || ir
->eI
== eiBD
)
463 fprintf(stderr
, "NOTE: The simulation uses pressure coupling and/or stochastic dynamics.\n"
464 "gmx convert-tpr can not provide binary identical continuation.\n"
465 "If you want that, supply a checkpoint file to mdrun\n\n");
468 if (EI_SD(ir
->eI
) || ir
->eI
== eiBD
)
470 fprintf(stderr
, "\nChanging ld-seed from %"GMX_PRId64
" ", ir
->ld_seed
);
471 ir
->ld_seed
= (gmx_int64_t
)gmx_rng_make_seed();
472 fprintf(stderr
, "to %"GMX_PRId64
"\n\n", ir
->ld_seed
);
475 frame_fn
= ftp2fn(efTRN
, NFILE
, fnm
);
477 if (fn2ftp(frame_fn
) == efCPT
)
482 "\nREADING STATE FROM CHECKPOINT %s...\n\n",
485 read_checkpoint_state(frame_fn
, &sim_part
,
486 &run_step
, &run_t
, &state
);
491 "\nREADING COORDS, VELS AND BOX FROM TRAJECTORY %s...\n\n",
494 fp
= open_trn(frame_fn
, "r");
497 fp_ener
= open_enx(ftp2fn(efEDR
, NFILE
, fnm
), "r");
498 do_enxnms(fp_ener
, &nre
, &enm
);
503 /* Now scan until the last set of x and v (step == 0)
504 * or the ones at step step.
510 bFrame
= fread_trnheader(fp
, &head
, &bOK
);
511 if (bOK
&& frame
== 0)
513 if (mtop
.natoms
!= head
.natoms
)
515 gmx_fatal(FARGS
, "Number of atoms in Topology (%d) "
516 "is not the same as in Trajectory (%d)\n",
517 mtop
.natoms
, head
.natoms
);
519 snew(newx
, head
.natoms
);
520 snew(newv
, head
.natoms
);
522 bFrame
= bFrame
&& bOK
;
525 bOK
= fread_htrn(fp
, &head
, newbox
, newx
, newv
, NULL
);
527 bFrame
= bFrame
&& bOK
;
530 (head
.x_size
) && (head
.v_size
|| !bVel
))
535 /* Read until the energy time is >= the trajectory time */
536 while (fr_ener
->t
< head
.t
&& do_enx(fp_ener
, fr_ener
))
540 bUse
= (fr_ener
->t
== head
.t
);
551 run_step
= head
.step
;
552 state
.fep_state
= head
.fep_state
;
553 state
.lambda
[efptFEP
] = head
.lambda
;
554 copy_mat(newbox
, state
.box
);
559 sprintf(buf
, "\r%s %s frame %s%s: step %s%s time %s",
560 "%s", "%s", "%6", GMX_PRId64
, "%6", GMX_PRId64
, " %8.3f");
562 bUse
? "Read " : "Skipped", ftp2ext(fn2ftp(frame_fn
)),
563 frame
, head
.step
, head
.t
);
565 if (bTime
&& (head
.t
>= start_t
))
574 free_enxframe(fr_ener
);
575 free_enxnms(nre
, enm
);
578 fprintf(stderr
, "\n");
582 fprintf(stderr
, "%s frame %s (step %s, time %g) is incomplete\n",
583 ftp2ext(fn2ftp(frame_fn
)), gmx_step_str(frame
-1, buf2
),
584 gmx_step_str(head
.step
, buf
), head
.t
);
586 fprintf(stderr
, "\nUsing frame of step %s time %g\n",
587 gmx_step_str(run_step
, buf
), run_t
);
593 get_enx_state(ftp2fn(efEDR
, NFILE
, fnm
), run_t
, &mtop
.groups
, ir
, &state
);
597 fprintf(stderr
, "\nWARNING: The simulation uses %s temperature and/or %s pressure coupling,\n"
598 " the continuation will only be exact when an energy file is supplied\n\n",
599 ETCOUPLTYPE(etcNOSEHOOVER
),
600 EPCOUPLTYPE(epcPARRINELLORAHMAN
));
605 ir
->fepvals
->init_fep_state
= init_fep_state
;
612 fprintf(stderr
, "Setting nsteps to %s\n", gmx_step_str(nsteps_req
, buf
));
613 ir
->nsteps
= nsteps_req
;
617 /* Determine total number of steps remaining */
620 ir
->nsteps
= ir
->nsteps
- (run_step
- ir
->init_step
) + (gmx_int64_t
)(extend_t
/ir
->delta_t
+ 0.5);
621 printf("Extending remaining runtime of by %g ps (now %s steps)\n",
622 extend_t
, gmx_step_str(ir
->nsteps
, buf
));
626 printf("nsteps = %s, run_step = %s, current_t = %g, until = %g\n",
627 gmx_step_str(ir
->nsteps
, buf
),
628 gmx_step_str(run_step
, buf2
),
630 ir
->nsteps
= (gmx_int64_t
)((until_t
- run_t
)/ir
->delta_t
+ 0.5);
631 printf("Extending remaining runtime until %g ps (now %s steps)\n",
632 until_t
, gmx_step_str(ir
->nsteps
, buf
));
636 ir
->nsteps
-= run_step
- ir
->init_step
;
638 printf("%s steps (%g ps) remaining from first run.\n",
639 gmx_step_str(ir
->nsteps
, buf
), ir
->nsteps
*ir
->delta_t
);
643 if (bNsteps
|| bZeroQ
|| (ir
->nsteps
> 0))
645 ir
->init_step
= run_step
;
647 if (ftp2bSet(efNDX
, NFILE
, fnm
) ||
648 !(bNsteps
|| bExtend
|| bUntil
|| bTraj
))
650 atoms
= gmx_mtop_global_atoms(&mtop
);
651 get_index(&atoms
, ftp2fn_null(efNDX
, NFILE
, fnm
), 1,
652 &gnx
, &index
, &grpname
);
655 bSel
= (gnx
!= state
.natoms
);
656 for (i
= 0; ((i
< gnx
) && (!bSel
)); i
++)
658 bSel
= (i
!= index
[i
]);
667 fprintf(stderr
, "Will write subset %s of original tpx containing %d "
668 "atoms\n", grpname
, gnx
);
669 reduce_topology_x(gnx
, index
, &mtop
, state
.x
, state
.v
);
675 fprintf(stderr
, "Zero-ing charges for group %s\n", grpname
);
679 fprintf(stderr
, "Will write full tpx file (no selection)\n");
683 state_t
= ir
->init_t
+ ir
->init_step
*ir
->delta_t
;
684 sprintf(buf
, "Writing statusfile with starting step %s%s and length %s%s steps...\n", "%10", GMX_PRId64
, "%10", GMX_PRId64
);
685 fprintf(stderr
, buf
, ir
->init_step
, ir
->nsteps
);
686 fprintf(stderr
, " time %10.3f and length %10.3f ps\n",
687 state_t
, ir
->nsteps
*ir
->delta_t
);
688 write_tpx_state(opt2fn("-o", NFILE
, fnm
), ir
, &state
, &mtop
);
692 printf("You've simulated long enough. Not writing tpr file\n");