Split espresso I/O routines from confio.cpp
[alexxy/gromacs.git] / src / gromacs / fileio / confio.cpp
index b9844ba9671bec3b1fd5a29134099ae3db9e189b..ac6390b2a79847f6c1ed2d66cdef4dad199f06f6 100644 (file)
 
 #include "confio.h"
 
-#include <cerrno>
-#include <cmath>
 #include <cstdio>
-#include <cstdlib>
 #include <cstring>
 
-#include <algorithm>
-
+#include "gromacs/fileio/espio.h"
 #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"
 #include "gromacs/topology/atoms.h"
 #include "gromacs/topology/block.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/smalloc.h"
 
-static int get_espresso_word(FILE *fp, char word[])
-{
-    int  ret, nc, i;
-
-    ret = 0;
-    nc  = 0;
-
-    do
-    {
-        i = fgetc(fp);
-        if (i != EOF)
-        {
-            if (i == ' ' || i == '\n' || i == '\t')
-            {
-                if (nc > 0)
-                {
-                    ret = 1;
-                }
-            }
-            else if (i == '{')
-            {
-                if (nc == 0)
-                {
-                    word[nc++] = '{';
-                }
-                ret = 2;
-            }
-            else if (i == '}')
-            {
-                if (nc == 0)
-                {
-                    word[nc++] = '}';
-                }
-                ret = 3;
-            }
-            else
-            {
-                word[nc++] = (char)i;
-            }
-        }
-    }
-    while (i != EOF && ret == 0);
-
-    word[nc] = '\0';
-
-    /*  printf("'%s'\n",word); */
-
-    return ret;
-}
-
-static int check_open_parenthesis(FILE *fp, int r,
-                                  const char *infile, const char *keyword)
-{
-    int  level_inc;
-    char word[STRLEN];
-
-    level_inc = 0;
-    if (r == 2)
-    {
-        level_inc++;
-    }
-    else
-    {
-        r = get_espresso_word(fp, word);
-        if (r == 2)
-        {
-            level_inc++;
-        }
-        else
-        {
-            gmx_fatal(FARGS, "Expected '{' after '%s' in file '%s'",
-                      keyword, infile);
-        }
-    }
-
-    return level_inc;
-}
-
-static int check_close_parenthesis(FILE *fp, int r,
-                                   const char *infile, const char *keyword)
-{
-    int  level_inc;
-    char word[STRLEN];
-
-    level_inc = 0;
-    if (r == 3)
-    {
-        level_inc--;
-    }
-    else
-    {
-        r = get_espresso_word(fp, word);
-        if (r == 3)
-        {
-            level_inc--;
-        }
-        else
-        {
-            gmx_fatal(FARGS, "Expected '}' after section '%s' in file '%s'",
-                      keyword, infile);
-        }
-    }
-
-    return level_inc;
-}
-
-enum {
-    espID, espPOS, espTYPE, espQ, espV, espF, espMOLECULE, espNR
-};
-const char *esp_prop[espNR] = {
-    "id", "pos", "type", "q", "v", "f",
-    "molecule"
-};
-
-static void read_espresso_conf(const char *infile, char *title,
-                               t_atoms *atoms, rvec x[], rvec *v, matrix box)
-{
-    t_symtab *symtab = NULL;
-    FILE     *fp;
-    char      word[STRLEN], buf[STRLEN];
-    int       level, r, nprop, p, i, m, molnr;
-    int       prop[32];
-    double    d;
-    gmx_bool  bFoundParticles, bFoundProp, bFoundVariable, bMol;
-
-    if (!symtab)
-    {
-        snew(symtab, 1);
-        open_symtab(symtab);
-    }
-    // TODO: The code does not understand titles it writes...
-    title[0] = '\0';
-
-    clear_mat(box);
-
-    fp = gmx_fio_fopen(infile, "r");
-
-    bFoundParticles = FALSE;
-    bFoundVariable  = FALSE;
-    bMol            = FALSE;
-    level           = 0;
-    while ((r = get_espresso_word(fp, word)))
-    {
-        if (level == 1 && std::strcmp(word, "particles") == 0 && !bFoundParticles)
-        {
-            bFoundParticles = TRUE;
-            level          += check_open_parenthesis(fp, r, infile, "particles");
-            nprop           = 0;
-            while (level == 2 && (r = get_espresso_word(fp, word)))
-            {
-                bFoundProp = FALSE;
-                for (p = 0; p < espNR; p++)
-                {
-                    if (strcmp(word, esp_prop[p]) == 0)
-                    {
-                        bFoundProp    = TRUE;
-                        prop[nprop++] = p;
-                        /* printf("  prop[%d] = %s\n",nprop-1,esp_prop[prop[nprop-1]]); */
-                    }
-                }
-                if (!bFoundProp && word[0] != '}')
-                {
-                    gmx_fatal(FARGS, "Can not read Espresso files with particle property '%s'", word);
-                }
-                if (bFoundProp && p == espMOLECULE)
-                {
-                    bMol = TRUE;
-                }
-                if (r == 3)
-                {
-                    level--;
-                }
-            }
-
-            i = 0;
-            while (level > 0 && (r = get_espresso_word(fp, word)))
-            {
-                if (r == 2)
-                {
-                    level++;
-                }
-                else if (r == 3)
-                {
-                    level--;
-                }
-                if (level == 2)
-                {
-                    for (p = 0; p < nprop; p++)
-                    {
-                        switch (prop[p])
-                        {
-                            case espID:
-                                r = get_espresso_word(fp, word);
-                                /* Not used */
-                                break;
-                            case espPOS:
-                                for (m = 0; m < 3; m++)
-                                {
-                                    r = get_espresso_word(fp, word);
-                                    sscanf(word, "%lf", &d);
-                                    x[i][m] = d;
-                                }
-                                break;
-                            case espTYPE:
-                                r                   = get_espresso_word(fp, word);
-                                atoms->atom[i].type = std::strtol(word, NULL, 10);
-                                break;
-                            case espQ:
-                                r = get_espresso_word(fp, word);
-                                sscanf(word, "%lf", &d);
-                                atoms->atom[i].q = d;
-                                break;
-                            case espV:
-                                for (m = 0; m < 3; m++)
-                                {
-                                    r = get_espresso_word(fp, word);
-                                    sscanf(word, "%lf", &d);
-                                    v[i][m] = d;
-                                }
-                                break;
-                            case espF:
-                                for (m = 0; m < 3; m++)
-                                {
-                                    r = get_espresso_word(fp, word);
-                                    /* not used */
-                                }
-                                break;
-                            case espMOLECULE:
-                                r     = get_espresso_word(fp, word);
-                                molnr = std::strtol(word, NULL, 10);
-                                if (i == 0 ||
-                                    atoms->resinfo[atoms->atom[i-1].resind].nr != molnr)
-                                {
-                                    atoms->atom[i].resind =
-                                        (i == 0 ? 0 : atoms->atom[i-1].resind+1);
-                                    atoms->resinfo[atoms->atom[i].resind].nr       = molnr;
-                                    atoms->resinfo[atoms->atom[i].resind].ic       = ' ';
-                                    atoms->resinfo[atoms->atom[i].resind].chainid  = ' ';
-                                    atoms->resinfo[atoms->atom[i].resind].chainnum = molnr; /* Not sure if this is right? */
-                                }
-                                else
-                                {
-                                    atoms->atom[i].resind = atoms->atom[i-1].resind;
-                                }
-                                break;
-                        }
-                    }
-                    /* Generate an atom name from the particle type */
-                    sprintf(buf, "T%d", atoms->atom[i].type);
-                    atoms->atomname[i] = put_symtab(symtab, buf);
-                    if (bMol)
-                    {
-                        if (i == 0 || atoms->atom[i].resind != atoms->atom[i-1].resind)
-                        {
-                            atoms->resinfo[atoms->atom[i].resind].name =
-                                put_symtab(symtab, "MOL");
-                        }
-                    }
-                    else
-                    {
-                        /* Residue number is the atom number */
-                        atoms->atom[i].resind = i;
-                        /* Generate an residue name from the particle type */
-                        if (atoms->atom[i].type < 26)
-                        {
-                            sprintf(buf, "T%c", 'A'+atoms->atom[i].type);
-                        }
-                        else
-                        {
-                            sprintf(buf, "T%c%c",
-                                    'A'+atoms->atom[i].type/26, 'A'+atoms->atom[i].type%26);
-                        }
-                        t_atoms_set_resinfo(atoms, i, symtab, buf, i, ' ', 0, ' ');
-                    }
-
-                    if (r == 3)
-                    {
-                        level--;
-                    }
-                    i++;
-                }
-            }
-            atoms->nres = atoms->nr;
-
-            if (i != atoms->nr)
-            {
-                gmx_fatal(FARGS, "Internal inconsistency in Espresso routines, read %d atoms, expected %d atoms", i, atoms->nr);
-            }
-        }
-        else if (level == 1 && std::strcmp(word, "variable") == 0 && !bFoundVariable)
-        {
-            bFoundVariable = TRUE;
-            level         += check_open_parenthesis(fp, r, infile, "variable");
-            while (level == 2 && (r = get_espresso_word(fp, word)))
-            {
-                if (level == 2 && std::strcmp(word, "box_l") == 0)
-                {
-                    for (m = 0; m < 3; m++)
-                    {
-                        r = get_espresso_word(fp, word);
-                        sscanf(word, "%lf", &d);
-                        box[m][m] = d;
-                    }
-                    level += check_close_parenthesis(fp, r, infile, "box_l");
-                }
-            }
-        }
-        else if (r == 2)
-        {
-            level++;
-        }
-        else if (r == 3)
-        {
-            level--;
-        }
-    }
-
-    if (!bFoundParticles)
-    {
-        fprintf(stderr, "Did not find a particles section in Espresso file '%s'\n",
-                infile);
-    }
-
-    gmx_fio_fclose(fp);
-}
-
-static int get_espresso_coordnum(const char *infile)
-{
-    FILE    *fp;
-    char     word[STRLEN];
-    int      natoms, level, r;
-    gmx_bool bFoundParticles;
-
-    natoms = 0;
-
-    fp = gmx_fio_fopen(infile, "r");
-
-    bFoundParticles = FALSE;
-    level           = 0;
-    while ((r = get_espresso_word(fp, word)) && !bFoundParticles)
-    {
-        if (level == 1 && strcmp(word, "particles") == 0 && !bFoundParticles)
-        {
-            bFoundParticles = TRUE;
-            level          += check_open_parenthesis(fp, r, infile, "particles");
-            while (level > 0 && (r = get_espresso_word(fp, word)))
-            {
-                if (r == 2)
-                {
-                    level++;
-                    if (level == 2)
-                    {
-                        natoms++;
-                    }
-                }
-                else if (r == 3)
-                {
-                    level--;
-                }
-            }
-        }
-        else if (r == 2)
-        {
-            level++;
-        }
-        else if (r == 3)
-        {
-            level--;
-        }
-    }
-    if (!bFoundParticles)
-    {
-        fprintf(stderr, "Did not find a particles section in Espresso file '%s'\n",
-                infile);
-    }
-
-    gmx_fio_fclose(fp);
-
-    return natoms;
-}
-
-static void write_espresso_conf_indexed(FILE *out, const char *title,
-                                        t_atoms *atoms, int nx, atom_id *index,
-                                        rvec *x, rvec *v, matrix box)
-{
-    int i, j;
-
-    fprintf(out, "# %s\n", title);
-    if (TRICLINIC(box))
-    {
-        gmx_warning("The Espresso format does not support triclinic unit-cells");
-    }
-    fprintf(out, "{variable {box_l %f %f %f}}\n", box[0][0], box[1][1], box[2][2]);
-
-    fprintf(out, "{particles {id pos type q%s}\n", v ? " v" : "");
-    for (i = 0; i < nx; i++)
-    {
-        if (index)
-        {
-            j = index[i];
-        }
-        else
-        {
-            j = i;
-        }
-        fprintf(out, "\t{%d %f %f %f %d %g",
-                j, x[j][XX], x[j][YY], x[j][ZZ],
-                atoms->atom[j].type, atoms->atom[j].q);
-        if (v)
-        {
-            fprintf(out, " %f %f %f", v[j][XX], v[j][YY], v[j][ZZ]);
-        }
-        fprintf(out, "}\n");
-    }
-    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 +130,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);