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,2015, 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/commandline/pargs.h"
42 #include "gromacs/fileio/enxio.h"
43 #include "gromacs/fileio/tpxio.h"
44 #include "gromacs/fileio/trrio.h"
45 #include "gromacs/gmxpreprocess/readir.h"
46 #include "gromacs/legacyheaders/checkpoint.h"
47 #include "gromacs/legacyheaders/names.h"
48 #include "gromacs/legacyheaders/types/ifunc.h"
49 #include "gromacs/legacyheaders/types/inputrec.h"
50 #include "gromacs/legacyheaders/types/simple.h"
51 #include "gromacs/legacyheaders/types/state.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/random/random.h"
54 #include "gromacs/topology/index.h"
55 #include "gromacs/topology/mtop_util.h"
56 #include "gromacs/topology/topology.h"
57 #include "gromacs/utility/arraysize.h"
58 #include "gromacs/utility/cstringutil.h"
59 #include "gromacs/utility/fatalerror.h"
60 #include "gromacs/utility/futil.h"
61 #include "gromacs/utility/smalloc.h"
63 #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))
65 static gmx_bool
*bKeepIt(int gnx
, int natoms
, atom_id index
[])
71 for (i
= 0; (i
< gnx
); i
++)
73 RANGECHK(index
[i
], natoms
);
80 static atom_id
*invind(int gnx
, int natoms
, atom_id index
[])
86 for (i
= 0; (i
< gnx
); i
++)
88 RANGECHK(index
[i
], natoms
);
95 static void reduce_block(gmx_bool bKeep
[], t_block
*block
,
101 snew(index
, block
->nr
);
104 for (i
= 0; (i
< block
->nr
); i
++)
106 for (j
= block
->index
[i
]; (j
< block
->index
[i
+1]); j
++)
113 if (newj
> index
[newi
])
120 fprintf(stderr
, "Reduced block %8s from %6d to %6d index-, %6d to %6d a-entries\n",
121 name
, block
->nr
, newi
, block
->index
[block
->nr
], newj
);
122 block
->index
= index
;
126 static void reduce_blocka(atom_id invindex
[], gmx_bool bKeep
[], t_blocka
*block
,
130 int i
, j
, k
, newi
, newj
;
132 snew(index
, block
->nr
);
136 for (i
= 0; (i
< block
->nr
); i
++)
138 for (j
= block
->index
[i
]; (j
< block
->index
[i
+1]); j
++)
143 a
[newj
] = invindex
[k
];
147 if (newj
> index
[newi
])
154 fprintf(stderr
, "Reduced block %8s from %6d to %6d index-, %6d to %6d a-entries\n",
155 name
, block
->nr
, newi
, block
->nra
, newj
);
156 block
->index
= index
;
162 static void reduce_rvec(int gnx
, atom_id index
[], rvec vv
[])
168 for (i
= 0; (i
< gnx
); i
++)
170 copy_rvec(vv
[index
[i
]], ptr
[i
]);
172 for (i
= 0; (i
< gnx
); i
++)
174 copy_rvec(ptr
[i
], vv
[i
]);
179 static void reduce_atom(int gnx
, atom_id index
[], t_atom atom
[], char ***atomname
,
180 int *nres
, t_resinfo
*resinfo
)
189 snew(rinfo
, atom
[index
[gnx
-1]].resind
+1);
190 for (i
= 0; (i
< gnx
); i
++)
192 ptr
[i
] = atom
[index
[i
]];
193 aname
[i
] = atomname
[index
[i
]];
196 for (i
= 0; (i
< gnx
); i
++)
199 atomname
[i
] = aname
[i
];
200 if ((i
== 0) || (atom
[i
].resind
!= atom
[i
-1].resind
))
203 rinfo
[nr
] = resinfo
[atom
[i
].resind
];
208 for (i
= 0; (i
< nr
); i
++)
210 resinfo
[i
] = rinfo
[i
];
219 static void reduce_ilist(atom_id invindex
[], gmx_bool bKeep
[],
220 t_ilist
*il
, int nratoms
, const char *name
)
230 for (i
= 0; (i
< il
->nr
); i
+= nratoms
+1)
233 for (j
= 1; (j
<= nratoms
); j
++)
235 bB
= bB
&& bKeep
[il
->iatoms
[i
+j
]];
239 ia
[newnr
++] = il
->iatoms
[i
];
240 for (j
= 1; (j
<= nratoms
); j
++)
242 ia
[newnr
++] = invindex
[il
->iatoms
[i
+j
]];
246 fprintf(stderr
, "Reduced ilist %8s from %6d to %6d entries\n",
247 name
, il
->nr
/(nratoms
+1),
251 for (i
= 0; (i
< newnr
); i
++)
253 il
->iatoms
[i
] = ia
[i
];
260 static void reduce_topology_x(int gnx
, atom_id index
[],
261 gmx_mtop_t
*mtop
, rvec x
[], rvec v
[])
268 top
= gmx_mtop_t_to_t_topology(mtop
);
269 bKeep
= bKeepIt(gnx
, top
.atoms
.nr
, index
);
270 invindex
= invind(gnx
, top
.atoms
.nr
, index
);
272 reduce_block(bKeep
, &(top
.cgs
), "cgs");
273 reduce_block(bKeep
, &(top
.mols
), "mols");
274 reduce_blocka(invindex
, bKeep
, &(top
.excls
), "excls");
275 reduce_rvec(gnx
, index
, x
);
276 reduce_rvec(gnx
, index
, v
);
277 reduce_atom(gnx
, index
, top
.atoms
.atom
, top
.atoms
.atomname
,
278 &(top
.atoms
.nres
), top
.atoms
.resinfo
);
280 for (i
= 0; (i
< F_NRE
); i
++)
282 reduce_ilist(invindex
, bKeep
, &(top
.idef
.il
[i
]),
283 interaction_function
[i
].nratoms
,
284 interaction_function
[i
].name
);
290 snew(mtop
->moltype
, mtop
->nmoltype
);
291 mtop
->moltype
[0].name
= mtop
->name
;
292 mtop
->moltype
[0].atoms
= top
.atoms
;
293 for (i
= 0; i
< F_NRE
; i
++)
295 mtop
->moltype
[0].ilist
[i
] = top
.idef
.il
[i
];
297 mtop
->moltype
[0].atoms
= top
.atoms
;
298 mtop
->moltype
[0].cgs
= top
.cgs
;
299 mtop
->moltype
[0].excls
= top
.excls
;
302 snew(mtop
->molblock
, mtop
->nmolblock
);
303 mtop
->molblock
[0].type
= 0;
304 mtop
->molblock
[0].nmol
= 1;
305 mtop
->molblock
[0].natoms_mol
= top
.atoms
.nr
;
306 mtop
->molblock
[0].nposres_xA
= 0;
307 mtop
->molblock
[0].nposres_xB
= 0;
309 mtop
->natoms
= top
.atoms
.nr
;
312 static void zeroq(atom_id index
[], gmx_mtop_t
*mtop
)
316 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
318 for (i
= 0; (i
< mtop
->moltype
[mt
].atoms
.nr
); i
++)
320 mtop
->moltype
[mt
].atoms
.atom
[index
[i
]].q
= 0;
321 mtop
->moltype
[mt
].atoms
.atom
[index
[i
]].qB
= 0;
326 int gmx_convert_tpr(int argc
, char *argv
[])
328 const char *desc
[] = {
329 "[THISMODULE] can edit run input files in four ways.[PAR]",
330 "[BB]1.[bb] by modifying the number of steps in a run input file",
331 "with options [TT]-extend[tt], [TT]-until[tt] or [TT]-nsteps[tt]",
332 "(nsteps=-1 means unlimited number of steps)[PAR]",
333 "[BB]2.[bb] (OBSOLETE) by creating a run input file",
334 "for a continuation run when your simulation has crashed due to e.g.",
335 "a full disk, or by making a continuation run input file.",
336 "This option is obsolete, since mdrun now writes and reads",
338 "[BB]Note[bb] that a frame with coordinates and velocities is needed.",
339 "When pressure and/or Nose-Hoover temperature coupling is used",
340 "an energy file can be supplied to get an exact continuation",
341 "of the original run.[PAR]",
342 "[BB]3.[bb] by creating a [REF].tpx[ref] file for a subset of your original",
343 "tpx file, which is useful when you want to remove the solvent from",
344 "your [REF].tpx[ref] file, or when you want to make e.g. a pure C[GRK]alpha[grk] [REF].tpx[ref] file.",
345 "Note that you may need to use [TT]-nsteps -1[tt] (or similar) to get",
347 "[BB]WARNING: this [REF].tpx[ref] file is not fully functional[bb].[PAR]",
348 "[BB]4.[bb] by setting the charges of a specified group",
349 "to zero. This is useful when doing free energy estimates",
350 "using the LIE (Linear Interaction Energy) method."
353 const char *top_fn
, *frame_fn
;
355 ener_file_t fp_ener
= NULL
;
356 gmx_trr_header_t head
;
358 gmx_int64_t nsteps_req
, run_step
, frame
;
359 double run_t
, state_t
;
360 gmx_bool bOK
, bNsteps
, bExtend
, bUntil
, bTime
, bTraj
;
361 gmx_bool bFrame
, bUse
, bSel
, bNeedEner
, bReadEner
, bScanEner
, bFepState
;
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 { efTPR
, NULL
, NULL
, ffREAD
},
378 { efTRN
, "-f", NULL
, ffOPTRD
},
379 { efEDR
, "-e", NULL
, ffOPTRD
},
380 { efNDX
, NULL
, NULL
, ffOPTRD
},
381 { efTPR
, "-o", "tprout", 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" },
409 /* Parse the command line */
410 if (!parse_common_args(&argc
, argv
, 0, NFILE
, fnm
, asize(pa
), pa
,
411 asize(desc
), desc
, 0, NULL
, &oenv
))
416 /* Convert int to gmx_int64_t */
417 nsteps_req
= nsteps_req_int
;
418 bNsteps
= opt2parg_bSet("-nsteps", asize(pa
), pa
);
419 bExtend
= opt2parg_bSet("-extend", asize(pa
), pa
);
420 bUntil
= opt2parg_bSet("-until", asize(pa
), pa
);
421 bFepState
= opt2parg_bSet("-init_fep_state", asize(pa
), pa
);
422 bTime
= opt2parg_bSet("-time", asize(pa
), pa
);
423 bTraj
= (opt2bSet("-f", NFILE
, fnm
) || bTime
);
425 top_fn
= ftp2fn(efTPR
, NFILE
, fnm
);
426 fprintf(stderr
, "Reading toplogy and stuff from %s\n", top_fn
);
429 read_tpx_state(top_fn
, ir
, &state
, &mtop
);
430 run_step
= ir
->init_step
;
431 run_t
= ir
->init_step
*ir
->delta_t
+ ir
->init_t
;
433 if (!EI_STATE_VELOCITY(ir
->eI
))
441 "NOTE: Reading the state from trajectory is an obsolete feature of gmx convert-tpr.\n"
442 " Continuation should be done by loading a checkpoint file with mdrun -cpi\n"
443 " This guarantees that all state variables are transferred.\n"
444 " gmx convert-tpr is now only useful for increasing nsteps,\n"
445 " but even that can often be avoided by using mdrun -maxh\n"
448 if (ir
->bContinuation
!= bContinuation
)
450 fprintf(stderr
, "Modifying ir->bContinuation to %s\n",
451 bool_names
[bContinuation
]);
453 ir
->bContinuation
= bContinuation
;
456 bNeedEner
= (ir
->epc
== epcPARRINELLORAHMAN
|| ir
->etc
== etcNOSEHOOVER
);
457 bReadEner
= (bNeedEner
&& ftp2bSet(efEDR
, NFILE
, fnm
));
458 bScanEner
= (bReadEner
&& !bTime
);
460 if (ir
->epc
!= epcNO
|| EI_SD(ir
->eI
) || ir
->eI
== eiBD
)
462 fprintf(stderr
, "NOTE: The simulation uses pressure coupling and/or stochastic dynamics.\n"
463 "gmx convert-tpr can not provide binary identical continuation.\n"
464 "If you want that, supply a checkpoint file to mdrun\n\n");
467 if (EI_SD(ir
->eI
) || ir
->eI
== eiBD
)
469 fprintf(stderr
, "\nChanging ld-seed from %" GMX_PRId64
" ", ir
->ld_seed
);
470 ir
->ld_seed
= (gmx_int64_t
)gmx_rng_make_seed();
471 fprintf(stderr
, "to %" GMX_PRId64
"\n\n", ir
->ld_seed
);
474 frame_fn
= ftp2fn(efTRN
, NFILE
, fnm
);
476 if (fn2ftp(frame_fn
) == efCPT
)
481 "\nREADING STATE FROM CHECKPOINT %s...\n\n",
484 read_checkpoint_state(frame_fn
, &sim_part
,
485 &run_step
, &run_t
, &state
);
490 "\nREADING COORDS, VELS AND BOX FROM TRAJECTORY %s...\n\n",
493 fp
= gmx_trr_open(frame_fn
, "r");
496 fp_ener
= open_enx(ftp2fn(efEDR
, NFILE
, fnm
), "r");
497 do_enxnms(fp_ener
, &nre
, &enm
);
502 /* Now scan until the last set of x and v (step == 0)
503 * or the ones at step step.
509 bFrame
= gmx_trr_read_frame_header(fp
, &head
, &bOK
);
510 if (bOK
&& frame
== 0)
512 if (mtop
.natoms
!= head
.natoms
)
514 gmx_fatal(FARGS
, "Number of atoms in Topology (%d) "
515 "is not the same as in Trajectory (%d)\n",
516 mtop
.natoms
, head
.natoms
);
518 snew(newx
, head
.natoms
);
519 snew(newv
, head
.natoms
);
521 bFrame
= bFrame
&& bOK
;
524 bOK
= gmx_trr_read_frame_data(fp
, &head
, newbox
, newx
, newv
, NULL
);
526 bFrame
= bFrame
&& bOK
;
529 (head
.x_size
) && (head
.v_size
|| !bVel
))
534 /* Read until the energy time is >= the trajectory time */
535 while (fr_ener
->t
< head
.t
&& do_enx(fp_ener
, fr_ener
))
539 bUse
= (fr_ener
->t
== head
.t
);
550 run_step
= head
.step
;
551 state
.fep_state
= head
.fep_state
;
552 state
.lambda
[efptFEP
] = head
.lambda
;
553 copy_mat(newbox
, state
.box
);
558 sprintf(buf
, "\r%s %s frame %s%s: step %s%s time %s",
559 "%s", "%s", "%6", GMX_PRId64
, "%6", GMX_PRId64
, " %8.3f");
561 bUse
? "Read " : "Skipped", ftp2ext(fn2ftp(frame_fn
)),
562 frame
, head
.step
, head
.t
);
564 if (bTime
&& (head
.t
>= start_t
))
573 free_enxframe(fr_ener
);
574 free_enxnms(nre
, enm
);
577 fprintf(stderr
, "\n");
581 fprintf(stderr
, "%s frame %s (step %s, time %g) is incomplete\n",
582 ftp2ext(fn2ftp(frame_fn
)), gmx_step_str(frame
-1, buf2
),
583 gmx_step_str(head
.step
, buf
), head
.t
);
585 fprintf(stderr
, "\nUsing frame of step %s time %g\n",
586 gmx_step_str(run_step
, buf
), run_t
);
592 get_enx_state(ftp2fn(efEDR
, NFILE
, fnm
), run_t
, &mtop
.groups
, ir
, &state
);
596 fprintf(stderr
, "\nWARNING: The simulation uses %s temperature and/or %s pressure coupling,\n"
597 " the continuation will only be exact when an energy file is supplied\n\n",
598 ETCOUPLTYPE(etcNOSEHOOVER
),
599 EPCOUPLTYPE(epcPARRINELLORAHMAN
));
604 ir
->fepvals
->init_fep_state
= init_fep_state
;
611 fprintf(stderr
, "Setting nsteps to %s\n", gmx_step_str(nsteps_req
, buf
));
612 ir
->nsteps
= nsteps_req
;
616 /* Determine total number of steps remaining */
619 ir
->nsteps
= ir
->nsteps
- (run_step
- ir
->init_step
) + (gmx_int64_t
)(extend_t
/ir
->delta_t
+ 0.5);
620 printf("Extending remaining runtime of by %g ps (now %s steps)\n",
621 extend_t
, gmx_step_str(ir
->nsteps
, buf
));
625 printf("nsteps = %s, run_step = %s, current_t = %g, until = %g\n",
626 gmx_step_str(ir
->nsteps
, buf
),
627 gmx_step_str(run_step
, buf2
),
629 ir
->nsteps
= (gmx_int64_t
)((until_t
- run_t
)/ir
->delta_t
+ 0.5);
630 printf("Extending remaining runtime until %g ps (now %s steps)\n",
631 until_t
, gmx_step_str(ir
->nsteps
, buf
));
635 ir
->nsteps
-= run_step
- ir
->init_step
;
637 printf("%s steps (%g ps) remaining from first run.\n",
638 gmx_step_str(ir
->nsteps
, buf
), ir
->nsteps
*ir
->delta_t
);
642 if (bNsteps
|| bZeroQ
|| (ir
->nsteps
> 0))
644 ir
->init_step
= run_step
;
646 if (ftp2bSet(efNDX
, NFILE
, fnm
) ||
647 !(bNsteps
|| bExtend
|| bUntil
|| bTraj
))
649 atoms
= gmx_mtop_global_atoms(&mtop
);
650 get_index(&atoms
, ftp2fn_null(efNDX
, NFILE
, fnm
), 1,
651 &gnx
, &index
, &grpname
);
654 bSel
= (gnx
!= state
.natoms
);
655 for (i
= 0; ((i
< gnx
) && (!bSel
)); i
++)
657 bSel
= (i
!= index
[i
]);
666 fprintf(stderr
, "Will write subset %s of original tpx containing %d "
667 "atoms\n", grpname
, gnx
);
668 reduce_topology_x(gnx
, index
, &mtop
, state
.x
, state
.v
);
674 fprintf(stderr
, "Zero-ing charges for group %s\n", grpname
);
678 fprintf(stderr
, "Will write full tpx file (no selection)\n");
682 state_t
= ir
->init_t
+ ir
->init_step
*ir
->delta_t
;
683 sprintf(buf
, "Writing statusfile with starting step %s%s and length %s%s steps...\n", "%10", GMX_PRId64
, "%10", GMX_PRId64
);
684 fprintf(stderr
, buf
, ir
->init_step
, ir
->nsteps
);
685 fprintf(stderr
, " time %10.3f and length %10.3f ps\n",
686 state_t
, ir
->nsteps
*ir
->delta_t
);
687 write_tpx_state(opt2fn("-o", NFILE
, fnm
), ir
, &state
, &mtop
);
691 printf("You've simulated long enough. Not writing tpr file\n");