Commit 2002cfde authored by Erik Lindahl's avatar Erik Lindahl Committed by David van der Spoel
Browse files

Improve frame time/step handling in trjconv

Store the exact step in PDB/GRO file headers,
and be more careful about not claiming to have
time or step information when it was not available.
This change will avoid some of the problems
described in #2189, but it does not yet properly
fix the issue in the TNG library.

Refs #2189.

Change-Id: I44bb59fbca83da53f6e8d4e494ae6476a82eb7cd
parent a726027a
......@@ -3,7 +3,7 @@
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
* Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
* Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
......@@ -385,6 +385,13 @@ gmx_bool gro_next_x_or_v(FILE *status, t_trxframe *fr)
}
}
if ((p = std::strstr(title, "step=")) != nullptr)
{
p += 5;
fr->step = 0; // Default value if fr-bStep is false
fr->bStep = (sscanf(p, "%" GMX_SCNd64, &fr->step) == 1);
}
if (atoms.nr != fr->natoms)
{
gmx_fatal(FARGS, "Number of atoms in gro frame (%d) doesn't match the number in the previous frame (%d)", atoms.nr, fr->natoms);
......
......@@ -1485,9 +1485,10 @@ gmx_bool gmx_read_next_tng_frame(tng_trajectory_t input,
fr->step = frameNumber;
fr->bStep = TRUE;
// Convert the time to ps
fr->time = frameTime / PICO;
fr->bTime = TRUE;
fr->bTime = (frameTime > 0);
// TODO This does not leak, but is not exception safe.
/* values must be freed before leaving this function */
......
......@@ -3,7 +3,7 @@
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
* Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
* Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
......@@ -725,7 +725,7 @@ static gmx_bool pdb_next_x(t_trxstatus *status, FILE *fp, t_trxframe *fr)
// It is not worthwhile introducing extra variables in the
// read_pdbfile call to verify that a model_nr was read.
int ePBC, model_nr = -1, na;
char title[STRLEN], *time;
char title[STRLEN], *time, *step;
double dbl;
atoms.nr = fr->natoms;
......@@ -751,32 +751,15 @@ static gmx_bool pdb_next_x(t_trxstatus *status, FILE *fp, t_trxframe *fr)
copy_mat(boxpdb, fr->box);
}
if (model_nr != -1)
{
fr->bStep = TRUE;
fr->step = model_nr;
}
time = std::strstr(title, " t= ");
if (time)
{
fr->bTime = TRUE;
sscanf(time+4, "%lf", &dbl);
fr->time = (real)dbl;
}
else
{
fr->bTime = FALSE;
/* this is a bit dirty, but it will work: if no time is read from
comment line in pdb file, set time to current frame number */
if (fr->bStep)
{
fr->time = (real)fr->step;
}
else
{
fr->time = (real)nframes_read(status);
}
}
fr->step = 0;
step = std::strstr(title, " step= ");
fr->bStep = (step && sscanf(step+7, "%" GMX_SCNd64, &fr->step) == 1);
dbl = 0.0;
time = std::strstr(title, " t= ");
fr->bTime = (time && sscanf(time+4, "%lf", &dbl) == 1);
fr->time = dbl;
if (na == 0)
{
return FALSE;
......
......@@ -3,7 +3,7 @@
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2007, The GROMACS development team.
* Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
* Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
......@@ -91,7 +91,9 @@ static void clust_size(const char *ndx, const char *trx, const char *xpm,
real tf, dx2, cut2, *t_x = nullptr, *t_y, cmid, cmax, cav, ekin;
int i, j, k, ai, aj, ci, cj, nframe, nclust, n_x, max_size = 0;
int *clust_index, *clust_size, max_clust_size, max_clust_ind, nav, nhisto;
t_rgb rlo = { 1.0, 1.0, 1.0 };
t_rgb rlo = { 1.0, 1.0, 1.0 };
int frameCounter = 0;
real frameTime;
clear_trxframe(&fr, TRUE);
auto timeLabel = output_env_get_time_label(oenv);
......@@ -265,7 +267,19 @@ static void clust_size(const char *ndx, const char *trx, const char *xpm,
}
n_x++;
srenew(t_x, n_x);
t_x[n_x-1] = fr.time*tf;
if (fr.bTime)
{
frameTime = fr.time;
}
else if (fr.bStep)
{
frameTime = fr.step;
}
else
{
frameTime = ++frameCounter;
}
t_x[n_x-1] = frameTime*tf;
srenew(cs_dist, n_x);
snew(cs_dist[n_x-1], nindex);
nclust = 0;
......@@ -291,12 +305,12 @@ static void clust_size(const char *ndx, const char *trx, const char *xpm,
}
}
}
fprintf(fp, "%14.6e %10d\n", fr.time, nclust);
fprintf(fp, "%14.6e %10d\n", frameTime, nclust);
if (nav > 0)
{
fprintf(gp, "%14.6e %10.3f\n", fr.time, cav/nav);
fprintf(gp, "%14.6e %10.3f\n", frameTime, cav/nav);
}
fprintf(hp, "%14.6e %10d\n", fr.time, max_clust_size);
fprintf(hp, "%14.6e %10d\n", frameTime, max_clust_size);
}
/* Analyse velocities, if present */
if (fr.bV)
......@@ -326,7 +340,7 @@ static void clust_size(const char *ndx, const char *trx, const char *xpm,
}
}
temp = (ekin*2.0)/(3.0*tfac*max_clust_size*BOLTZ);
fprintf(tp, "%10.3f %10.3f\n", fr.time, temp);
fprintf(tp, "%10.3f %10.3f\n", frameTime, temp);
}
}
}
......
......@@ -904,7 +904,7 @@ int gmx_trjconv(int argc, char *argv[])
char out_file2[256], *charpt;
char *outf_base = nullptr;
const char *outf_ext = nullptr;
char top_title[256], title[256], filemode[5];
char top_title[256], title[256], timestr[32], stepstr[32], filemode[5];
gmx_output_env_t *oenv;
t_filenm fnm[] = {
......@@ -1105,13 +1105,20 @@ int gmx_trjconv(int argc, char *argv[])
}
/* top_title is only used for gro and pdb,
* the header in such a file is top_title t= ...
* to prevent a double t=, remove it from top_title
* the header in such a file is top_title, followed by
* t= ... and/or step= ...
* to prevent double t= or step=, remove it from top_title.
* From GROMACS-2018 we only write t/step when the frame actually
* has a valid time/step, so we need to check for both separately.
*/
if ((charpt = std::strstr(top_title, " t= ")))
{
charpt[0] = '\0';
}
if ((charpt = std::strstr(top_title, " step= ")))
{
charpt[0] = '\0';
}
if (bCONECT)
{
......@@ -1827,8 +1834,29 @@ int gmx_trjconv(int argc, char *argv[])
case efGRO:
case efG96:
case efPDB:
sprintf(title, "Generated by trjconv : %s t= %9.5f",
top_title, frout.time);
// Only add a generator statement if title is empty,
// to avoid multiple generated-by statements from various programs
if (std::strlen(top_title) == 0)
{
sprintf(top_title, "Generated by trjconv");
}
if (frout.bTime)
{
sprintf(timestr, " t= %9.5f", frout.time);
}
else
{
std::strcpy(timestr, "");
}
if (frout.bStep)
{
sprintf(stepstr, " step= %" GMX_PRId64, frout.step);
}
else
{
std::strcpy(stepstr, "");
}
snprintf(title, 256, "%s%s%s", top_title, timestr, stepstr);
if (bSeparate || bSplitHere)
{
out = gmx_ffopen(out_file2, "w");
......@@ -1841,7 +1869,6 @@ int gmx_trjconv(int argc, char *argv[])
break;
case efPDB:
fprintf(out, "REMARK GENERATED BY TRJCONV\n");
sprintf(title, "%s t= %9.5f", top_title, frout.time);
/* if reading from pdb, we want to keep the original
model numbering else we write the output frame
number plus one, because model 0 is not allowed in pdb */
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment