Split gro I/O routines from confio.cpp
authorTeemu Murtola <teemu.murtola@gmail.com>
Wed, 22 Jul 2015 18:08:39 +0000 (21:08 +0300)
committerTeemu Murtola <teemu.murtola@gmail.com>
Wed, 22 Jul 2015 19:02:16 +0000 (21:02 +0200)
Move routines to do gro-format I/O to a separate file.  As with g96io.h,
this would ideally be internal to the fileio module, but unfortunately a
few tools peek into the details of some specific format.

The naming of the routines will be changed separately.

Change-Id: Ia0b190a95f4fb2cd223a6581531aaf45e99d3440

src/gromacs/fileio/confio.cpp
src/gromacs/fileio/confio.h
src/gromacs/fileio/groio.cpp [new file with mode: 0644]
src/gromacs/fileio/groio.h [new file with mode: 0644]
src/gromacs/fileio/trxio.cpp
src/gromacs/gmxana/gmx_confrms.c
src/gromacs/gmxana/gmx_trjconv.c

index b9844ba9671bec3b1fd5a29134099ae3db9e189b..14093492279b5a84534198799be21ea7e873b7e3 100644 (file)
 #include "gromacs/fileio/filenm.h"
 #include "gromacs/fileio/g96io.h"
 #include "gromacs/fileio/gmxfio.h"
+#include "gromacs/fileio/groio.h"
 #include "gromacs/fileio/pdbio.h"
 #include "gromacs/fileio/tpxio.h"
 #include "gromacs/fileio/trx.h"
 #include "gromacs/fileio/trxio.h"
-#include "gromacs/legacyheaders/copyrite.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/topology/atomprop.h"
@@ -483,549 +483,6 @@ static void write_espresso_conf_indexed(FILE *out, const char *title,
     fprintf(out, "}\n");
 }
 
-static void get_coordnum_fp (FILE *in, char *title, int *natoms)
-{
-    char line[STRLEN+1];
-
-    fgets2 (title, STRLEN, in);
-    fgets2 (line, STRLEN, in);
-    if (sscanf (line, "%d", natoms) != 1)
-    {
-        gmx_fatal(FARGS, "gro file does not have the number of atoms on the second line");
-    }
-}
-
-static void get_coordnum (const char *infile, int *natoms)
-{
-    FILE *in;
-    char  title[STRLEN];
-
-    in = gmx_fio_fopen(infile, "r");
-    get_coordnum_fp(in, title, natoms);
-    gmx_fio_fclose (in);
-}
-
-static gmx_bool get_w_conf(FILE *in, const char *infile, char *title,
-                           t_symtab *symtab, t_atoms *atoms, int *ndec,
-                           rvec x[], rvec *v, matrix box)
-{
-    char       name[6];
-    char       resname[6], oldresname[6];
-    char       line[STRLEN+1], *ptr;
-    char       buf[256];
-    double     x1, y1, z1, x2, y2, z2;
-    rvec       xmin, xmax;
-    int        natoms, i, m, resnr, newres, oldres, ddist, c;
-    gmx_bool   bFirst, bVel;
-    char      *p1, *p2, *p3;
-
-    newres  = -1;
-    oldres  = -12345; /* Unlikely number for the first residue! */
-    ddist   = 0;
-
-    /* Read the title and number of atoms */
-    get_coordnum_fp(in, title, &natoms);
-
-    if (natoms > atoms->nr)
-    {
-        gmx_fatal(FARGS, "gro file contains more atoms (%d) than expected (%d)",
-                  natoms, atoms->nr);
-    }
-    else if (natoms <  atoms->nr)
-    {
-        fprintf(stderr, "Warning: gro file contains less atoms (%d) than expected"
-                " (%d)\n", natoms, atoms->nr);
-    }
-
-    bFirst = TRUE;
-
-    bVel = FALSE;
-
-    resname[0]     = '\0';
-    oldresname[0]  = '\0';
-
-    /* just pray the arrays are big enough */
-    for (i = 0; (i < natoms); i++)
-    {
-        if ((fgets2 (line, STRLEN, in)) == NULL)
-        {
-            gmx_fatal(FARGS, "Unexpected end of file in file %s at line %d",
-                      infile, i+2);
-        }
-        if (strlen(line) < 39)
-        {
-            gmx_fatal(FARGS, "Invalid line in %s for atom %d:\n%s", infile, i+1, line);
-        }
-
-        /* determine read precision from distance between periods
-           (decimal points) */
-        if (bFirst)
-        {
-            bFirst = FALSE;
-            p1     = strchr(line, '.');
-            if (p1 == NULL)
-            {
-                gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
-            }
-            p2 = strchr(&p1[1], '.');
-            if (p2 == NULL)
-            {
-                gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
-            }
-            ddist = p2 - p1;
-            *ndec = ddist - 5;
-
-            p3 = strchr(&p2[1], '.');
-            if (p3 == NULL)
-            {
-                gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
-            }
-
-            if (p3 - p2 != ddist)
-            {
-                gmx_fatal(FARGS, "The spacing of the decimal points in file %s is not consistent for x, y and z", infile);
-            }
-        }
-
-        /* residue number*/
-        memcpy(name, line, 5);
-        name[5] = '\0';
-        sscanf(name, "%d", &resnr);
-        sscanf(line+5, "%5s", resname);
-
-        if (resnr != oldres || strncmp(resname, oldresname, sizeof(resname)))
-        {
-            oldres = resnr;
-            newres++;
-            if (newres >= natoms)
-            {
-                gmx_fatal(FARGS, "More residues than atoms in %s (natoms = %d)",
-                          infile, natoms);
-            }
-            atoms->atom[i].resind = newres;
-            t_atoms_set_resinfo(atoms, i, symtab, resname, resnr, ' ', 0, ' ');
-        }
-        else
-        {
-            atoms->atom[i].resind = newres;
-        }
-
-        /* atomname */
-        std::memcpy(name, line+10, 5);
-        atoms->atomname[i] = put_symtab(symtab, name);
-
-        /* Copy resname to oldresname after we are done with the sanity check above */
-        std::strncpy(oldresname, resname, sizeof(oldresname));
-
-        /* eventueel controle atomnumber met i+1 */
-
-        /* coordinates (start after residue data) */
-        ptr = line + 20;
-        /* Read fixed format */
-        for (m = 0; m < DIM; m++)
-        {
-            for (c = 0; (c < ddist && ptr[0]); c++)
-            {
-                buf[c] = ptr[0];
-                ptr++;
-            }
-            buf[c] = '\0';
-            if (sscanf (buf, "%lf %lf", &x1, &x2) != 1)
-            {
-                gmx_fatal(FARGS, "Something is wrong in the coordinate formatting of file %s. Note that gro is fixed format (see the manual)", infile);
-            }
-            else
-            {
-                x[i][m] = x1;
-            }
-        }
-
-        /* velocities (start after residues and coordinates) */
-        if (v)
-        {
-            /* Read fixed format */
-            for (m = 0; m < DIM; m++)
-            {
-                for (c = 0; (c < ddist && ptr[0]); c++)
-                {
-                    buf[c] = ptr[0];
-                    ptr++;
-                }
-                buf[c] = '\0';
-                if (sscanf (buf, "%lf", &x1) != 1)
-                {
-                    v[i][m] = 0;
-                }
-                else
-                {
-                    v[i][m] = x1;
-                    bVel    = TRUE;
-                }
-            }
-        }
-    }
-    atoms->nres = newres + 1;
-
-    /* box */
-    fgets2 (line, STRLEN, in);
-    if (sscanf (line, "%lf%lf%lf", &x1, &y1, &z1) != 3)
-    {
-        gmx_warning("Bad box in file %s", infile);
-
-        /* Generate a cubic box */
-        for (m = 0; (m < DIM); m++)
-        {
-            xmin[m] = xmax[m] = x[0][m];
-        }
-        for (i = 1; (i < atoms->nr); i++)
-        {
-            for (m = 0; (m < DIM); m++)
-            {
-                xmin[m] = std::min(xmin[m], x[i][m]);
-                xmax[m] = std::max(xmax[m], x[i][m]);
-            }
-        }
-        for (i = 0; i < DIM; i++)
-        {
-            for (m = 0; m < DIM; m++)
-            {
-                box[i][m] = 0.0;
-            }
-        }
-        for (m = 0; (m < DIM); m++)
-        {
-            box[m][m] = (xmax[m]-xmin[m]);
-        }
-        fprintf(stderr, "Generated a cubic box %8.3f x %8.3f x %8.3f\n",
-                box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
-    }
-    else
-    {
-        /* We found the first three values, the diagonal elements */
-        box[XX][XX] = x1;
-        box[YY][YY] = y1;
-        box[ZZ][ZZ] = z1;
-        if (sscanf (line, "%*f%*f%*f%lf%lf%lf%lf%lf%lf",
-                    &x1, &y1, &z1, &x2, &y2, &z2) != 6)
-        {
-            x1 = y1 = z1 = x2 = y2 = z2 = 0.0;
-        }
-        box[XX][YY] = x1;
-        box[XX][ZZ] = y1;
-        box[YY][XX] = z1;
-        box[YY][ZZ] = x2;
-        box[ZZ][XX] = y2;
-        box[ZZ][YY] = z2;
-    }
-
-    return bVel;
-}
-
-static void read_whole_conf(const char *infile, char *title,
-                            t_atoms *atoms, rvec x[], rvec *v, matrix box)
-{
-    FILE    *in;
-    int      ndec;
-    t_symtab symtab;
-
-    /* open file */
-    in = gmx_fio_fopen(infile, "r");
-
-    open_symtab(&symtab);
-    get_w_conf(in, infile, title, &symtab, atoms, &ndec, x, v, box);
-    /* We can't free the symbols, as they are still used in atoms, so
-     * the only choice is to leak them. */
-    free_symtab(&symtab);
-
-    gmx_fio_fclose(in);
-}
-
-static gmx_bool gmx_one_before_eof(FILE *fp)
-{
-    char     data[4];
-    gmx_bool beof;
-
-    if ((beof = fread(data, 1, 1, fp)) == 1)
-    {
-        gmx_fseek(fp, -1, SEEK_CUR);
-    }
-    return !beof;
-}
-
-gmx_bool gro_next_x_or_v(FILE *status, t_trxframe *fr)
-{
-    t_atoms  atoms;
-    t_symtab symtab;
-    char     title[STRLEN], *p;
-    double   tt;
-    int      ndec = 0, i;
-
-    if (gmx_one_before_eof(status))
-    {
-        return FALSE;
-    }
-
-    open_symtab(&symtab);
-    atoms.nr = fr->natoms;
-    snew(atoms.atom, fr->natoms);
-    atoms.nres = fr->natoms;
-    snew(atoms.resinfo, fr->natoms);
-    snew(atoms.atomname, fr->natoms);
-
-    fr->bV    = get_w_conf(status, title, title, &symtab, &atoms, &ndec, fr->x, fr->v, fr->box);
-    fr->bPrec = TRUE;
-    fr->prec  = 1;
-    /* prec = 10^ndec: */
-    for (i = 0; i < ndec; i++)
-    {
-        fr->prec *= 10;
-    }
-    fr->title  = title;
-    fr->bTitle = TRUE;
-    fr->bX     = TRUE;
-    fr->bBox   = TRUE;
-
-    sfree(atoms.atom);
-    sfree(atoms.resinfo);
-    sfree(atoms.atomname);
-    done_symtab(&symtab);
-
-    if ((p = strstr(title, "t=")) != NULL)
-    {
-        p += 2;
-        if (sscanf(p, "%lf", &tt) == 1)
-        {
-            fr->time  = tt;
-            fr->bTime = TRUE;
-        }
-        else
-        {
-            fr->time  = 0;
-            fr->bTime = FALSE;
-        }
-    }
-
-    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);
-    }
-
-    return TRUE;
-}
-
-int gro_first_x_or_v(FILE *status, t_trxframe *fr)
-{
-    char title[STRLEN];
-
-    frewind(status);
-    fprintf(stderr, "Reading frames from gro file");
-    get_coordnum_fp(status, title, &fr->natoms);
-    frewind(status);
-    fprintf(stderr, " '%s', %d atoms.\n", title, fr->natoms);
-    fr->bTitle = TRUE;
-    fr->title  = title;
-    if (fr->natoms == 0)
-    {
-        gmx_file("No coordinates in gro file");
-    }
-
-    snew(fr->x, fr->natoms);
-    snew(fr->v, fr->natoms);
-    gro_next_x_or_v(status, fr);
-
-    return fr->natoms;
-}
-
-static void make_hconf_format(int pr, gmx_bool bVel, char format[])
-{
-    int l, vpr;
-
-    /* build format string for printing,
-       something like "%8.3f" for x and "%8.4f" for v */
-    if (pr < 0)
-    {
-        pr = 0;
-    }
-    if (pr > 30)
-    {
-        pr = 30;
-    }
-    l   = pr+5;
-    vpr = pr+1;
-    if (bVel)
-    {
-        sprintf(format, "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
-                l, pr, l, pr, l, pr, l, vpr, l, vpr, l, vpr);
-    }
-    else
-    {
-        sprintf(format, "%%%d.%df%%%d.%df%%%d.%df\n", l, pr, l, pr, l, pr);
-    }
-
-}
-
-static void write_hconf_box(FILE *out, int pr, matrix box)
-{
-    char format[100];
-    int  l;
-
-    if (pr < 5)
-    {
-        pr = 5;
-    }
-    l = pr+5;
-
-    if (box[XX][YY] || box[XX][ZZ] || box[YY][XX] || box[YY][ZZ] ||
-        box[ZZ][XX] || box[ZZ][YY])
-    {
-        sprintf(format, "%%%d.%df%%%d.%df%%%d.%df"
-                "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
-                l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr);
-        fprintf(out, format,
-                box[XX][XX], box[YY][YY], box[ZZ][ZZ],
-                box[XX][YY], box[XX][ZZ], box[YY][XX],
-                box[YY][ZZ], box[ZZ][XX], box[ZZ][YY]);
-    }
-    else
-    {
-        sprintf(format, "%%%d.%df%%%d.%df%%%d.%df\n", l, pr, l, pr, l, pr);
-        fprintf(out, format,
-                box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
-    }
-}
-
-void write_hconf_indexed_p(FILE *out, const char *title, t_atoms *atoms,
-                           int nx, const atom_id index[], int pr,
-                           rvec *x, rvec *v, matrix box)
-{
-    char resnm[6], nm[6], format[100];
-    int  ai, i, resind, resnr;
-
-    bromacs(format, 99);
-    fprintf (out, "%s\n", (title && title[0]) ? title : format);
-    fprintf (out, "%5d\n", nx);
-
-    make_hconf_format(pr, v != NULL, format);
-
-    for (i = 0; (i < nx); i++)
-    {
-        ai = index[i];
-
-        resind = atoms->atom[ai].resind;
-        std::strncpy(resnm, " ??? ", sizeof(resnm)-1);
-        if (resind < atoms->nres)
-        {
-            std::strncpy(resnm, *atoms->resinfo[resind].name, sizeof(resnm)-1);
-            resnr = atoms->resinfo[resind].nr;
-        }
-        else
-        {
-            std::strncpy(resnm, " ??? ", sizeof(resnm)-1);
-            resnr = resind + 1;
-        }
-
-        if (atoms->atom)
-        {
-            std::strncpy(nm, *atoms->atomname[ai], sizeof(nm)-1);
-        }
-        else
-        {
-            std::strncpy(nm, " ??? ", sizeof(nm)-1);
-        }
-
-        fprintf(out, "%5d%-5.5s%5.5s%5d", resnr%100000, resnm, nm, (ai+1)%100000);
-        /* next fprintf uses built format string */
-        if (v)
-        {
-            fprintf(out, format,
-                    x[ai][XX], x[ai][YY], x[ai][ZZ], v[ai][XX], v[ai][YY], v[ai][ZZ]);
-        }
-        else
-        {
-            fprintf(out, format,
-                    x[ai][XX], x[ai][YY], x[ai][ZZ]);
-        }
-    }
-
-    write_hconf_box(out, pr, box);
-
-    fflush(out);
-}
-
-static void write_hconf_mtop(FILE *out, const char *title, gmx_mtop_t *mtop,
-                             int pr,
-                             rvec *x, rvec *v, matrix box)
-{
-    char                    format[100];
-    int                     i, resnr;
-    gmx_mtop_atomloop_all_t aloop;
-    t_atom                 *atom;
-    char                   *atomname, *resname;
-
-    bromacs(format, 99);
-    fprintf (out, "%s\n", (title && title[0]) ? title : format);
-    fprintf (out, "%5d\n", mtop->natoms);
-
-    make_hconf_format(pr, v != NULL, format);
-
-    aloop = gmx_mtop_atomloop_all_init(mtop);
-    while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
-    {
-        gmx_mtop_atomloop_all_names(aloop, &atomname, &resnr, &resname);
-
-        fprintf(out, "%5d%-5.5s%5.5s%5d",
-                resnr%100000, resname, atomname, (i+1)%100000);
-        /* next fprintf uses built format string */
-        if (v)
-        {
-            fprintf(out, format,
-                    x[i][XX], x[i][YY], x[i][ZZ], v[i][XX], v[i][YY], v[i][ZZ]);
-        }
-        else
-        {
-            fprintf(out, format,
-                    x[i][XX], x[i][YY], x[i][ZZ]);
-        }
-    }
-
-    write_hconf_box(out, pr, box);
-
-    fflush(out);
-}
-
-void write_hconf_p(FILE *out, const char *title, t_atoms *atoms, int pr,
-                   rvec *x, rvec *v, matrix box)
-{
-    atom_id *aa;
-    int      i;
-
-    snew(aa, atoms->nr);
-    for (i = 0; (i < atoms->nr); i++)
-    {
-        aa[i] = i;
-    }
-    write_hconf_indexed_p(out, title, atoms, atoms->nr, aa, pr, x, v, box);
-    sfree(aa);
-}
-
-void write_conf_p(const char *outfile, const char *title,
-                  t_atoms *atoms, int pr,
-                  rvec *x, rvec *v, matrix box)
-{
-    FILE *out;
-
-    out = gmx_fio_fopen(outfile, "w");
-    write_hconf_p(out, title, atoms, pr, x, v, box);
-
-    gmx_fio_fclose (out);
-}
-
-static void write_conf(const char *outfile, const char *title, t_atoms *atoms,
-                       rvec *x, rvec *v, matrix box)
-{
-    write_conf_p(outfile, title, atoms, 3, x, v, box);
-}
-
 void write_sto_conf_indexed(const char *outfile, const char *title,
                             t_atoms *atoms,
                             rvec x[], rvec *v, int ePBC, matrix box,
@@ -1095,7 +552,7 @@ void write_sto_conf(const char *outfile, const char *title, t_atoms *atoms,
     switch (ftp)
     {
         case efGRO:
-            write_conf(outfile, title, atoms, x, v, box);
+            write_conf_p(outfile, title, atoms, 3, x, v, box);
             break;
         case efG96:
             clear_trxframe(&fr, TRUE);
index 194c514e515b85484139b4c9b887559ed4cfc496..2c5dc2dc0e8ff2c80033708c6d47f6ade0f3fa38 100644 (file)
@@ -39,7 +39,6 @@
 
 #include <stdio.h>
 
-#include "gromacs/fileio/trx.h"
 #include "gromacs/legacyheaders/types/simple.h"
 
 /* For reading coordinate files it is assumed that enough memory
@@ -53,19 +52,6 @@ struct gmx_mtop_t;
 struct t_atoms;
 struct t_topology;
 
-gmx_bool gro_next_x_or_v(FILE *status, t_trxframe *fr);
-int gro_first_x_or_v(FILE *status, t_trxframe *fr);
-/* read first/next x and/or v frame from gro file */
-
-void write_hconf_indexed_p(FILE *out, const char *title, struct t_atoms *atoms,
-                           int nx, const atom_id index[], int ndec,
-                           rvec *x, rvec *v, matrix box);
-
-void write_hconf_p(FILE *out, const char *title, struct t_atoms *atoms, int ndec,
-                   rvec *x, rvec *v, matrix box);
-/* Write a Gromos file with precision ndec: number of decimal places in x,
- * v has one place more. */
-
 void write_sto_conf_indexed(const char *outfile, const char *title,
                             struct t_atoms *atoms,
                             rvec x[], rvec *v, int ePBC, matrix box,
@@ -106,4 +92,4 @@ gmx_bool read_tps_conf(const char *infile, char *title, struct t_topology *top,
 }
 #endif
 
-#endif  /* GMX_FILEIO_CONFIO_H */
+#endif
diff --git a/src/gromacs/fileio/groio.cpp b/src/gromacs/fileio/groio.cpp
new file mode 100644 (file)
index 0000000..47e65dd
--- /dev/null
@@ -0,0 +1,591 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014,2015, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#include "gmxpre.h"
+
+#include "groio.h"
+
+#include <cstdio>
+#include <cstring>
+
+#include <algorithm>
+
+#include "gromacs/fileio/gmxfio.h"
+#include "gromacs/fileio/trx.h"
+#include "gromacs/legacyheaders/copyrite.h"
+#include "gromacs/topology/atoms.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/topology/symtab.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/smalloc.h"
+
+static void get_coordnum_fp(FILE *in, char *title, int *natoms)
+{
+    char line[STRLEN+1];
+
+    fgets2(title, STRLEN, in);
+    fgets2(line, STRLEN, in);
+    if (sscanf(line, "%d", natoms) != 1)
+    {
+        gmx_fatal(FARGS, "gro file does not have the number of atoms on the second line");
+    }
+}
+
+void get_coordnum(const char *infile, int *natoms)
+{
+    FILE *in;
+    char  title[STRLEN];
+
+    in = gmx_fio_fopen(infile, "r");
+    get_coordnum_fp(in, title, natoms);
+    gmx_fio_fclose(in);
+}
+
+static gmx_bool get_w_conf(FILE *in, const char *infile, char *title,
+                           t_symtab *symtab, t_atoms *atoms, int *ndec,
+                           rvec x[], rvec *v, matrix box)
+{
+    char       name[6];
+    char       resname[6], oldresname[6];
+    char       line[STRLEN+1], *ptr;
+    char       buf[256];
+    double     x1, y1, z1, x2, y2, z2;
+    rvec       xmin, xmax;
+    int        natoms, i, m, resnr, newres, oldres, ddist, c;
+    gmx_bool   bFirst, bVel;
+    char      *p1, *p2, *p3;
+
+    newres  = -1;
+    oldres  = -12345; /* Unlikely number for the first residue! */
+    ddist   = 0;
+
+    /* Read the title and number of atoms */
+    get_coordnum_fp(in, title, &natoms);
+
+    if (natoms > atoms->nr)
+    {
+        gmx_fatal(FARGS, "gro file contains more atoms (%d) than expected (%d)",
+                  natoms, atoms->nr);
+    }
+    else if (natoms <  atoms->nr)
+    {
+        fprintf(stderr, "Warning: gro file contains less atoms (%d) than expected"
+                " (%d)\n", natoms, atoms->nr);
+    }
+
+    bFirst = TRUE;
+
+    bVel = FALSE;
+
+    resname[0]     = '\0';
+    oldresname[0]  = '\0';
+
+    /* just pray the arrays are big enough */
+    for (i = 0; (i < natoms); i++)
+    {
+        if ((fgets2(line, STRLEN, in)) == NULL)
+        {
+            gmx_fatal(FARGS, "Unexpected end of file in file %s at line %d",
+                      infile, i+2);
+        }
+        if (strlen(line) < 39)
+        {
+            gmx_fatal(FARGS, "Invalid line in %s for atom %d:\n%s", infile, i+1, line);
+        }
+
+        /* determine read precision from distance between periods
+           (decimal points) */
+        if (bFirst)
+        {
+            bFirst = FALSE;
+            p1     = strchr(line, '.');
+            if (p1 == NULL)
+            {
+                gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
+            }
+            p2 = strchr(&p1[1], '.');
+            if (p2 == NULL)
+            {
+                gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
+            }
+            ddist = p2 - p1;
+            *ndec = ddist - 5;
+
+            p3 = strchr(&p2[1], '.');
+            if (p3 == NULL)
+            {
+                gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
+            }
+
+            if (p3 - p2 != ddist)
+            {
+                gmx_fatal(FARGS, "The spacing of the decimal points in file %s is not consistent for x, y and z", infile);
+            }
+        }
+
+        /* residue number*/
+        memcpy(name, line, 5);
+        name[5] = '\0';
+        sscanf(name, "%d", &resnr);
+        sscanf(line+5, "%5s", resname);
+
+        if (resnr != oldres || strncmp(resname, oldresname, sizeof(resname)))
+        {
+            oldres = resnr;
+            newres++;
+            if (newres >= natoms)
+            {
+                gmx_fatal(FARGS, "More residues than atoms in %s (natoms = %d)",
+                          infile, natoms);
+            }
+            atoms->atom[i].resind = newres;
+            t_atoms_set_resinfo(atoms, i, symtab, resname, resnr, ' ', 0, ' ');
+        }
+        else
+        {
+            atoms->atom[i].resind = newres;
+        }
+
+        /* atomname */
+        std::memcpy(name, line+10, 5);
+        atoms->atomname[i] = put_symtab(symtab, name);
+
+        /* Copy resname to oldresname after we are done with the sanity check above */
+        std::strncpy(oldresname, resname, sizeof(oldresname));
+
+        /* eventueel controle atomnumber met i+1 */
+
+        /* coordinates (start after residue data) */
+        ptr = line + 20;
+        /* Read fixed format */
+        for (m = 0; m < DIM; m++)
+        {
+            for (c = 0; (c < ddist && ptr[0]); c++)
+            {
+                buf[c] = ptr[0];
+                ptr++;
+            }
+            buf[c] = '\0';
+            if (sscanf(buf, "%lf %lf", &x1, &x2) != 1)
+            {
+                gmx_fatal(FARGS, "Something is wrong in the coordinate formatting of file %s. Note that gro is fixed format (see the manual)", infile);
+            }
+            else
+            {
+                x[i][m] = x1;
+            }
+        }
+
+        /* velocities (start after residues and coordinates) */
+        if (v)
+        {
+            /* Read fixed format */
+            for (m = 0; m < DIM; m++)
+            {
+                for (c = 0; (c < ddist && ptr[0]); c++)
+                {
+                    buf[c] = ptr[0];
+                    ptr++;
+                }
+                buf[c] = '\0';
+                if (sscanf(buf, "%lf", &x1) != 1)
+                {
+                    v[i][m] = 0;
+                }
+                else
+                {
+                    v[i][m] = x1;
+                    bVel    = TRUE;
+                }
+            }
+        }
+    }
+    atoms->nres = newres + 1;
+
+    /* box */
+    fgets2(line, STRLEN, in);
+    if (sscanf(line, "%lf%lf%lf", &x1, &y1, &z1) != 3)
+    {
+        gmx_warning("Bad box in file %s", infile);
+
+        /* Generate a cubic box */
+        for (m = 0; (m < DIM); m++)
+        {
+            xmin[m] = xmax[m] = x[0][m];
+        }
+        for (i = 1; (i < atoms->nr); i++)
+        {
+            for (m = 0; (m < DIM); m++)
+            {
+                xmin[m] = std::min(xmin[m], x[i][m]);
+                xmax[m] = std::max(xmax[m], x[i][m]);
+            }
+        }
+        for (i = 0; i < DIM; i++)
+        {
+            for (m = 0; m < DIM; m++)
+            {
+                box[i][m] = 0.0;
+            }
+        }
+        for (m = 0; (m < DIM); m++)
+        {
+            box[m][m] = (xmax[m]-xmin[m]);
+        }
+        fprintf(stderr, "Generated a cubic box %8.3f x %8.3f x %8.3f\n",
+                box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
+    }
+    else
+    {
+        /* We found the first three values, the diagonal elements */
+        box[XX][XX] = x1;
+        box[YY][YY] = y1;
+        box[ZZ][ZZ] = z1;
+        if (sscanf (line, "%*f%*f%*f%lf%lf%lf%lf%lf%lf",
+                    &x1, &y1, &z1, &x2, &y2, &z2) != 6)
+        {
+            x1 = y1 = z1 = x2 = y2 = z2 = 0.0;
+        }
+        box[XX][YY] = x1;
+        box[XX][ZZ] = y1;
+        box[YY][XX] = z1;
+        box[YY][ZZ] = x2;
+        box[ZZ][XX] = y2;
+        box[ZZ][YY] = z2;
+    }
+
+    return bVel;
+}
+
+void read_whole_conf(const char *infile, char *title,
+                     t_atoms *atoms, rvec x[], rvec *v, matrix box)
+{
+    FILE    *in;
+    int      ndec;
+    t_symtab symtab;
+
+    /* open file */
+    in = gmx_fio_fopen(infile, "r");
+
+    open_symtab(&symtab);
+    get_w_conf(in, infile, title, &symtab, atoms, &ndec, x, v, box);
+    /* We can't free the symbols, as they are still used in atoms, so
+     * the only choice is to leak them. */
+    free_symtab(&symtab);
+
+    gmx_fio_fclose(in);
+}
+
+static gmx_bool gmx_one_before_eof(FILE *fp)
+{
+    char     data[4];
+    gmx_bool beof;
+
+    if ((beof = fread(data, 1, 1, fp)) == 1)
+    {
+        gmx_fseek(fp, -1, SEEK_CUR);
+    }
+    return !beof;
+}
+
+gmx_bool gro_next_x_or_v(FILE *status, t_trxframe *fr)
+{
+    t_atoms  atoms;
+    t_symtab symtab;
+    char     title[STRLEN], *p;
+    double   tt;
+    int      ndec = 0, i;
+
+    if (gmx_one_before_eof(status))
+    {
+        return FALSE;
+    }
+
+    open_symtab(&symtab);
+    atoms.nr = fr->natoms;
+    snew(atoms.atom, fr->natoms);
+    atoms.nres = fr->natoms;
+    snew(atoms.resinfo, fr->natoms);
+    snew(atoms.atomname, fr->natoms);
+
+    fr->bV    = get_w_conf(status, title, title, &symtab, &atoms, &ndec, fr->x, fr->v, fr->box);
+    fr->bPrec = TRUE;
+    fr->prec  = 1;
+    /* prec = 10^ndec: */
+    for (i = 0; i < ndec; i++)
+    {
+        fr->prec *= 10;
+    }
+    fr->title  = title;
+    fr->bTitle = TRUE;
+    fr->bX     = TRUE;
+    fr->bBox   = TRUE;
+
+    sfree(atoms.atom);
+    sfree(atoms.resinfo);
+    sfree(atoms.atomname);
+    done_symtab(&symtab);
+
+    if ((p = strstr(title, "t=")) != NULL)
+    {
+        p += 2;
+        if (sscanf(p, "%lf", &tt) == 1)
+        {
+            fr->time  = tt;
+            fr->bTime = TRUE;
+        }
+        else
+        {
+            fr->time  = 0;
+            fr->bTime = FALSE;
+        }
+    }
+
+    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);
+    }
+
+    return TRUE;
+}
+
+int gro_first_x_or_v(FILE *status, t_trxframe *fr)
+{
+    char title[STRLEN];
+
+    frewind(status);
+    fprintf(stderr, "Reading frames from gro file");
+    get_coordnum_fp(status, title, &fr->natoms);
+    frewind(status);
+    fprintf(stderr, " '%s', %d atoms.\n", title, fr->natoms);
+    fr->bTitle = TRUE;
+    fr->title  = title;
+    if (fr->natoms == 0)
+    {
+        gmx_file("No coordinates in gro file");
+    }
+
+    snew(fr->x, fr->natoms);
+    snew(fr->v, fr->natoms);
+    gro_next_x_or_v(status, fr);
+
+    return fr->natoms;
+}
+
+static void make_hconf_format(int pr, gmx_bool bVel, char format[])
+{
+    int l, vpr;
+
+    /* build format string for printing,
+       something like "%8.3f" for x and "%8.4f" for v */
+    if (pr < 0)
+    {
+        pr = 0;
+    }
+    if (pr > 30)
+    {
+        pr = 30;
+    }
+    l   = pr+5;
+    vpr = pr+1;
+    if (bVel)
+    {
+        sprintf(format, "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
+                l, pr, l, pr, l, pr, l, vpr, l, vpr, l, vpr);
+    }
+    else
+    {
+        sprintf(format, "%%%d.%df%%%d.%df%%%d.%df\n", l, pr, l, pr, l, pr);
+    }
+
+}
+
+static void write_hconf_box(FILE *out, int pr, matrix box)
+{
+    char format[100];
+    int  l;
+
+    if (pr < 5)
+    {
+        pr = 5;
+    }
+    l = pr+5;
+
+    if (box[XX][YY] || box[XX][ZZ] || box[YY][XX] || box[YY][ZZ] ||
+        box[ZZ][XX] || box[ZZ][YY])
+    {
+        sprintf(format, "%%%d.%df%%%d.%df%%%d.%df"
+                "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
+                l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr);
+        fprintf(out, format,
+                box[XX][XX], box[YY][YY], box[ZZ][ZZ],
+                box[XX][YY], box[XX][ZZ], box[YY][XX],
+                box[YY][ZZ], box[ZZ][XX], box[ZZ][YY]);
+    }
+    else
+    {
+        sprintf(format, "%%%d.%df%%%d.%df%%%d.%df\n", l, pr, l, pr, l, pr);
+        fprintf(out, format,
+                box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
+    }
+}
+
+void write_hconf_indexed_p(FILE *out, const char *title, t_atoms *atoms,
+                           int nx, const atom_id index[], int pr,
+                           rvec *x, rvec *v, matrix box)
+{
+    char resnm[6], nm[6], format[100];
+    int  ai, i, resind, resnr;
+
+    bromacs(format, 99);
+    fprintf(out, "%s\n", (title && title[0]) ? title : format);
+    fprintf(out, "%5d\n", nx);
+
+    make_hconf_format(pr, v != NULL, format);
+
+    for (i = 0; (i < nx); i++)
+    {
+        ai = index[i];
+
+        resind = atoms->atom[ai].resind;
+        std::strncpy(resnm, " ??? ", sizeof(resnm)-1);
+        if (resind < atoms->nres)
+        {
+            std::strncpy(resnm, *atoms->resinfo[resind].name, sizeof(resnm)-1);
+            resnr = atoms->resinfo[resind].nr;
+        }
+        else
+        {
+            std::strncpy(resnm, " ??? ", sizeof(resnm)-1);
+            resnr = resind + 1;
+        }
+
+        if (atoms->atom)
+        {
+            std::strncpy(nm, *atoms->atomname[ai], sizeof(nm)-1);
+        }
+        else
+        {
+            std::strncpy(nm, " ??? ", sizeof(nm)-1);
+        }
+
+        fprintf(out, "%5d%-5.5s%5.5s%5d", resnr%100000, resnm, nm, (ai+1)%100000);
+        /* next fprintf uses built format string */
+        if (v)
+        {
+            fprintf(out, format,
+                    x[ai][XX], x[ai][YY], x[ai][ZZ], v[ai][XX], v[ai][YY], v[ai][ZZ]);
+        }
+        else
+        {
+            fprintf(out, format,
+                    x[ai][XX], x[ai][YY], x[ai][ZZ]);
+        }
+    }
+
+    write_hconf_box(out, pr, box);
+
+    fflush(out);
+}
+
+void write_hconf_mtop(FILE *out, const char *title, gmx_mtop_t *mtop, int pr,
+                      rvec *x, rvec *v, matrix box)
+{
+    char                    format[100];
+    int                     i, resnr;
+    gmx_mtop_atomloop_all_t aloop;
+    t_atom                 *atom;
+    char                   *atomname, *resname;
+
+    bromacs(format, 99);
+    fprintf(out, "%s\n", (title && title[0]) ? title : format);
+    fprintf(out, "%5d\n", mtop->natoms);
+
+    make_hconf_format(pr, v != NULL, format);
+
+    aloop = gmx_mtop_atomloop_all_init(mtop);
+    while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
+    {
+        gmx_mtop_atomloop_all_names(aloop, &atomname, &resnr, &resname);
+
+        fprintf(out, "%5d%-5.5s%5.5s%5d",
+                resnr%100000, resname, atomname, (i+1)%100000);
+        /* next fprintf uses built format string */
+        if (v)
+        {
+            fprintf(out, format,
+                    x[i][XX], x[i][YY], x[i][ZZ], v[i][XX], v[i][YY], v[i][ZZ]);
+        }
+        else
+        {
+            fprintf(out, format,
+                    x[i][XX], x[i][YY], x[i][ZZ]);
+        }
+    }
+
+    write_hconf_box(out, pr, box);
+
+    fflush(out);
+}
+
+void write_hconf_p(FILE *out, const char *title, t_atoms *atoms, int pr,
+                   rvec *x, rvec *v, matrix box)
+{
+    atom_id *aa;
+    int      i;
+
+    snew(aa, atoms->nr);
+    for (i = 0; (i < atoms->nr); i++)
+    {
+        aa[i] = i;
+    }
+    write_hconf_indexed_p(out, title, atoms, atoms->nr, aa, pr, x, v, box);
+    sfree(aa);
+}
+
+void write_conf_p(const char *outfile, const char *title,
+                  t_atoms *atoms, int pr,
+                  rvec *x, rvec *v, matrix box)
+{
+    FILE *out;
+
+    out = gmx_fio_fopen(outfile, "w");
+    write_hconf_p(out, title, atoms, pr, x, v, box);
+    gmx_fio_fclose(out);
+}
diff --git a/src/gromacs/fileio/groio.h b/src/gromacs/fileio/groio.h
new file mode 100644 (file)
index 0000000..b89c8ee
--- /dev/null
@@ -0,0 +1,81 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014,2015, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#ifndef GMX_FILEIO_GROIO_H
+#define GMX_FILEIO_GROIO_H
+
+#include <stdio.h>
+
+#include "gromacs/legacyheaders/types/simple.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct gmx_mtop_t;
+struct t_atoms;
+struct t_topology;
+struct t_trxframe;
+
+void get_coordnum(const char *infile, int *natoms);
+void read_whole_conf(const char *infile, char *title,
+                     struct t_atoms *atoms, rvec x[], rvec *v, matrix box);
+
+gmx_bool gro_next_x_or_v(FILE *status, struct t_trxframe *fr);
+int gro_first_x_or_v(FILE *status, struct t_trxframe *fr);
+/* read first/next x and/or v frame from gro file */
+
+void write_hconf_indexed_p(FILE *out, const char *title, struct t_atoms *atoms,
+                           int nx, const atom_id index[], int ndec,
+                           rvec *x, rvec *v, matrix box);
+
+void write_hconf_mtop(FILE *out, const char *title, struct gmx_mtop_t *mtop, int pr,
+                      rvec *x, rvec *v, matrix box);
+
+void write_hconf_p(FILE *out, const char *title, struct t_atoms *atoms, int ndec,
+                   rvec *x, rvec *v, matrix box);
+/* Write a Gromos file with precision ndec: number of decimal places in x,
+ * v has one place more. */
+
+void write_conf_p(const char *outfile, const char *title,
+                  struct t_atoms *atoms, int pr,
+                  rvec *x, rvec *v, matrix box);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
index addaed90d829404e5feee305eb1050dfd11bbb1c..116dc4c2a89a50633f9688ead9a19b15997d8d04 100644 (file)
@@ -47,6 +47,7 @@
 #include "gromacs/fileio/g96io.h"
 #include "gromacs/fileio/gmxfio.h"
 #include "gromacs/fileio/gmxfio-xdr.h"
+#include "gromacs/fileio/groio.h"
 #include "gromacs/fileio/pdbio.h"
 #include "gromacs/fileio/timecontrol.h"
 #include "gromacs/fileio/tngio.h"
index b4489a7e67ee7ca5060e22781f31ae6df80dcb1a..26add192a6f9ad1c628a5b2578557dfa24160305 100644 (file)
@@ -42,6 +42,7 @@
 #include "gromacs/commandline/pargs.h"
 #include "gromacs/fileio/confio.h"
 #include "gromacs/fileio/filenm.h"
+#include "gromacs/fileio/groio.h"
 #include "gromacs/fileio/pdbio.h"
 #include "gromacs/gmxana/gmx_ana.h"
 #include "gromacs/legacyheaders/macros.h"
index eb83b301d410bd91d1353c431064655169607a4a..9db84bee50493dd295a364001c5e80fbb7aaec1e 100644 (file)
@@ -44,6 +44,7 @@
 #include "gromacs/fileio/confio.h"
 #include "gromacs/fileio/g96io.h"
 #include "gromacs/fileio/gmxfio.h"
+#include "gromacs/fileio/groio.h"
 #include "gromacs/fileio/pdbio.h"
 #include "gromacs/fileio/tngio_for_tools.h"
 #include "gromacs/fileio/tpxio.h"