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

diff --git a/src/gromacs/fileio/groio.cpp b/src/gromacs/fileio/groio.cpp
index 08fa47ddc896032a884c7dac141862d93bb3e9fa..a616f2a0418f3665b283e40b7ea6992547f9918a 100644 (file)
--- a/src/gromacs/fileio/groio.cpp
+++ b/src/gromacs/fileio/groio.cpp
@@ -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);

diff --git a/src/gromacs/fileio/tngio.cpp b/src/gromacs/fileio/tngio.cpp
index db99f290692fbd2dfe275e6507a80eaa8e08c7e7..b26ec181ecb3ed633f101f419cac08a8c372be13 100644 (file)
--- a/src/gromacs/fileio/tngio.cpp
+++ b/src/gromacs/fileio/tngio.cpp
@@ -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 */

diff --git a/src/gromacs/fileio/trxio.cpp b/src/gromacs/fileio/trxio.cpp
index 20ce0c01739733294ea578086c28b82242d47ab1..2dd6c05b4fe7fcf4781c1fca75c87ab158149a0f 100644 (file)
--- a/src/gromacs/fileio/trxio.cpp
+++ b/src/gromacs/fileio/trxio.cpp
@@ -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;

diff --git a/src/gromacs/gmxana/gmx_clustsize.cpp b/src/gromacs/gmxana/gmx_clustsize.cpp
index a0d6d507934af14c8ca2949619b3745c9a460b65..f1eb7d508c065b06b4c748cdb0426fea8cdb4a32 100644 (file)
--- a/src/gromacs/gmxana/gmx_clustsize.cpp
+++ b/src/gromacs/gmxana/gmx_clustsize.cpp
@@ -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);
                 }
             }
         }

diff --git a/src/gromacs/gmxana/gmx_trjconv.cpp b/src/gromacs/gmxana/gmx_trjconv.cpp
index e49700f45439b81d39ef3d0884015fc85c9c0d09..fd6ed925759b211fb1d418c208f5d8fdccc00d4d 100644 (file)
--- a/src/gromacs/gmxana/gmx_trjconv.cpp
+++ b/src/gromacs/gmxana/gmx_trjconv.cpp
@@ -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 */