Split g96 I/O routines from confio.cpp
[alexxy/gromacs.git] / src / gromacs / fileio / confio.cpp
index a7d9fdaa5d483d4eaf70282fc4927d6c4cec8096..b9844ba9671bec3b1fd5a29134099ae3db9e189b 100644 (file)
 #include <algorithm>
 
 #include "gromacs/fileio/filenm.h"
+#include "gromacs/fileio/g96io.h"
 #include "gromacs/fileio/gmxfio.h"
 #include "gromacs/fileio/pdbio.h"
 #include "gromacs/fileio/tpxio.h"
 #include "gromacs/fileio/trx.h"
 #include "gromacs/fileio/trxio.h"
-#include "gromacs/fileio/xdrf.h"
 #include "gromacs/legacyheaders/copyrite.h"
-#include "gromacs/legacyheaders/macros.h"
-#include "gromacs/legacyheaders/typedefs.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/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"
 
-#define CHAR_SHIFT 24
-
-static int read_g96_pos(char line[], t_symtab *symtab,
-                        FILE *fp, const char *infile,
-                        t_trxframe *fr)
-{
-    t_atoms   *atoms;
-    gmx_bool   bEnd;
-    int        nwanted, natoms, atnr, resnr = 0, oldres, newres, shift;
-    char       anm[STRLEN], resnm[STRLEN];
-    char       c1, c2;
-    double     db1, db2, db3;
-
-    nwanted = fr->natoms;
-
-    atoms = fr->atoms;
-
-    natoms = 0;
-
-    if (fr->bX)
-    {
-        if (fr->bAtoms)
-        {
-            shift = CHAR_SHIFT;
-        }
-        else
-        {
-            shift = 0;
-        }
-        newres  = -1;
-        oldres  = -666; /* Unlikely number for the first residue! */
-        bEnd    = FALSE;
-        while (!bEnd && fgets2(line, STRLEN, fp))
-        {
-            bEnd = (std::strncmp(line, "END", 3) == 0);
-            if (!bEnd  && (line[0] != '#'))
-            {
-                if (sscanf(line+shift, "%15lf%15lf%15lf", &db1, &db2, &db3) != 3)
-                {
-                    gmx_fatal(FARGS, "Did not find 3 coordinates for atom %d in %s\n",
-                              natoms+1, infile);
-                }
-                if ((nwanted != -1) && (natoms >= nwanted))
-                {
-                    gmx_fatal(FARGS,
-                              "Found more coordinates (%d) in %s than expected %d\n",
-                              natoms, infile, nwanted);
-                }
-                if (atoms)
-                {
-                    if (fr->bAtoms &&
-                        (sscanf(line, "%5d%c%5s%c%5s%7d", &resnr, &c1, resnm, &c2, anm, &atnr)
-                         != 6))
-                    {
-                        if (oldres >= 0)
-                        {
-                            resnr = oldres;
-                        }
-                        else
-                        {
-                            resnr    = 1;
-                            strncpy(resnm, "???", sizeof(resnm)-1);
-                        }
-                        strncpy(anm, "???", sizeof(anm)-1);
-                    }
-                    atoms->atomname[natoms] = put_symtab(symtab, anm);
-                    if (resnr != oldres)
-                    {
-                        oldres = resnr;
-                        newres++;
-                        if (newres >= atoms->nr)
-                        {
-                            gmx_fatal(FARGS, "More residues than atoms in %s (natoms = %d)",
-                                      infile, atoms->nr);
-                        }
-                        atoms->atom[natoms].resind = newres;
-                        if (newres+1 > atoms->nres)
-                        {
-                            atoms->nres = newres+1;
-                        }
-                        t_atoms_set_resinfo(atoms, natoms, symtab, resnm, resnr, ' ', 0, ' ');
-                    }
-                    else
-                    {
-                        atoms->atom[natoms].resind = newres;
-                    }
-                }
-                if (fr->x)
-                {
-                    fr->x[natoms][0] = db1;
-                    fr->x[natoms][1] = db2;
-                    fr->x[natoms][2] = db3;
-                }
-                natoms++;
-            }
-        }
-        if ((nwanted != -1) && natoms != nwanted)
-        {
-            fprintf(stderr,
-                    "Warning: found less coordinates (%d) in %s than expected %d\n",
-                    natoms, infile, nwanted);
-        }
-    }
-
-    fr->natoms = natoms;
-
-    return natoms;
-}
-
-static int read_g96_vel(char line[], FILE *fp, const char *infile,
-                        t_trxframe *fr)
-{
-    gmx_bool   bEnd;
-    int        nwanted, natoms = -1, shift;
-    double     db1, db2, db3;
-
-    nwanted = fr->natoms;
-
-    if (fr->v && fr->bV)
-    {
-        if (strcmp(line, "VELOCITYRED") == 0)
-        {
-            shift = 0;
-        }
-        else
-        {
-            shift = CHAR_SHIFT;
-        }
-        natoms = 0;
-        bEnd   = FALSE;
-        while (!bEnd && fgets2(line, STRLEN, fp))
-        {
-            bEnd = (strncmp(line, "END", 3) == 0);
-            if (!bEnd && (line[0] != '#'))
-            {
-                if (sscanf(line+shift, "%15lf%15lf%15lf", &db1, &db2, &db3) != 3)
-                {
-                    gmx_fatal(FARGS, "Did not find 3 velocities for atom %d in %s",
-                              natoms+1, infile);
-                }
-                if ((nwanted != -1) && (natoms >= nwanted))
-                {
-                    gmx_fatal(FARGS, "Found more velocities (%d) in %s than expected %d\n",
-                              natoms, infile, nwanted);
-                }
-                if (fr->v)
-                {
-                    fr->v[natoms][0] = db1;
-                    fr->v[natoms][1] = db2;
-                    fr->v[natoms][2] = db3;
-                }
-                natoms++;
-            }
-        }
-        if ((nwanted != -1) && (natoms != nwanted))
-        {
-            fprintf(stderr,
-                    "Warning: found less velocities (%d) in %s than expected %d\n",
-                    natoms, infile, nwanted);
-        }
-    }
-
-    return natoms;
-}
-
-int read_g96_conf(FILE *fp, const char *infile, t_trxframe *fr, char *line)
-{
-    t_symtab  *symtab = NULL;
-    gmx_bool   bAtStart, bTime, bAtoms, bPos, bVel, bBox, bEnd, bFinished;
-    int        natoms, nbp;
-    double     db1, db2, db3, db4, db5, db6, db7, db8, db9;
-
-    bAtStart = (ftell(fp) == 0);
-
-    clear_trxframe(fr, FALSE);
-
-    if (!symtab)
-    {
-        snew(symtab, 1);
-        open_symtab(symtab);
-    }
-
-    natoms = 0;
-
-    if (bAtStart)
-    {
-        while (!fr->bTitle && fgets2(line, STRLEN, fp))
-        {
-            fr->bTitle = (std::strcmp(line, "TITLE") == 0);
-        }
-        if (fr->title == NULL)
-        {
-            fgets2(line, STRLEN, fp);
-            fr->title = gmx_strdup(line);
-        }
-        bEnd = FALSE;
-        while (!bEnd && fgets2(line, STRLEN, fp))
-        {
-            bEnd = (std::strcmp(line, "END") == 0);
-        }
-        fgets2(line, STRLEN, fp);
-    }
-
-    /* Do not get a line if we are not at the start of the file, *
-     * because without a parameter file we don't know what is in *
-     * the trajectory and we have already read the line in the   *
-     * previous call (VERY DIRTY).                               */
-    bFinished = FALSE;
-    do
-    {
-        bTime  = (std::strcmp(line, "TIMESTEP") == 0);
-        bAtoms = (std::strcmp(line, "POSITION") == 0);
-        bPos   = (bAtoms || (strcmp(line, "POSITIONRED") == 0));
-        bVel   = (std::strncmp(line, "VELOCITY", 8) == 0);
-        bBox   = (std::strcmp(line, "BOX") == 0);
-        if (bTime)
-        {
-            if (!fr->bTime && !fr->bX)
-            {
-                fr->bStep = bTime;
-                fr->bTime = bTime;
-                do
-                {
-                    bFinished = (fgets2(line, STRLEN, fp) == NULL);
-                }
-                while (!bFinished && (line[0] == '#'));
-                sscanf(line, "%15d%15lf", &(fr->step), &db1);
-                fr->time = db1;
-            }
-            else
-            {
-                bFinished = TRUE;
-            }
-        }
-        if (bPos)
-        {
-            if (!fr->bX)
-            {
-                fr->bAtoms = bAtoms;
-                fr->bX     = bPos;
-                natoms     = read_g96_pos(line, symtab, fp, infile, fr);
-            }
-            else
-            {
-                bFinished = TRUE;
-            }
-        }
-        if (fr->v && bVel)
-        {
-            fr->bV = bVel;
-            natoms = read_g96_vel(line, fp, infile, fr);
-        }
-        if (bBox)
-        {
-            fr->bBox = bBox;
-            clear_mat(fr->box);
-            bEnd = FALSE;
-            while (!bEnd && fgets2(line, STRLEN, fp))
-            {
-                bEnd = (strncmp(line, "END", 3) == 0);
-                if (!bEnd && (line[0] != '#'))
-                {
-                    nbp = sscanf(line, "%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf",
-                                 &db1, &db2, &db3, &db4, &db5, &db6, &db7, &db8, &db9);
-                    if (nbp < 3)
-                    {
-                        gmx_fatal(FARGS, "Found a BOX line, but no box in %s", infile);
-                    }
-                    fr->box[XX][XX] = db1;
-                    fr->box[YY][YY] = db2;
-                    fr->box[ZZ][ZZ] = db3;
-                    if (nbp == 9)
-                    {
-                        fr->box[XX][YY] = db4;
-                        fr->box[XX][ZZ] = db5;
-                        fr->box[YY][XX] = db6;
-                        fr->box[YY][ZZ] = db7;
-                        fr->box[ZZ][XX] = db8;
-                        fr->box[ZZ][YY] = db9;
-                    }
-                }
-            }
-            bFinished = TRUE;
-        }
-    }
-    while (!bFinished && fgets2(line, STRLEN, fp));
-
-    free_symtab(symtab);
-
-    fr->natoms = natoms;
-
-    return natoms;
-}
-
-void write_g96_conf(FILE *out, t_trxframe *fr,
-                    int nindex, const atom_id *index)
-{
-    t_atoms *atoms;
-    int      nout, i, a;
-
-    atoms = fr->atoms;
-
-    if (index)
-    {
-        nout = nindex;
-    }
-    else
-    {
-        nout = fr->natoms;
-    }
-
-    if (fr->bTitle)
-    {
-        fprintf(out, "TITLE\n%s\nEND\n", fr->title);
-    }
-    if (fr->bStep || fr->bTime)
-    {
-        /* Officially the time format is %15.9, which is not enough for 10 ns */
-        fprintf(out, "TIMESTEP\n%15d%15.6f\nEND\n", fr->step, fr->time);
-    }
-    if (fr->bX)
-    {
-        if (fr->bAtoms)
-        {
-            fprintf(out, "POSITION\n");
-            for (i = 0; i < nout; i++)
-            {
-                if (index)
-                {
-                    a = index[i];
-                }
-                else
-                {
-                    a = i;
-                }
-                fprintf(out, "%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
-                        (atoms->resinfo[atoms->atom[a].resind].nr) % 100000,
-                        *atoms->resinfo[atoms->atom[a].resind].name,
-                        *atoms->atomname[a], (i+1) % 10000000,
-                        fr->x[a][XX], fr->x[a][YY], fr->x[a][ZZ]);
-            }
-        }
-        else
-        {
-            fprintf(out, "POSITIONRED\n");
-            for (i = 0; i < nout; i++)
-            {
-                if (index)
-                {
-                    a = index[i];
-                }
-                else
-                {
-                    a = i;
-                }
-                fprintf(out, "%15.9f%15.9f%15.9f\n",
-                        fr->x[a][XX], fr->x[a][YY], fr->x[a][ZZ]);
-            }
-        }
-        fprintf(out, "END\n");
-    }
-    if (fr->bV)
-    {
-        if (fr->bAtoms)
-        {
-            fprintf(out, "VELOCITY\n");
-            for (i = 0; i < nout; i++)
-            {
-                if (index)
-                {
-                    a = index[i];
-                }
-                else
-                {
-                    a = i;
-                }
-                fprintf(out, "%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
-                        (atoms->resinfo[atoms->atom[a].resind].nr) % 100000,
-                        *atoms->resinfo[atoms->atom[a].resind].name,
-                        *atoms->atomname[a], (i+1) % 10000000,
-                        fr->v[a][XX], fr->v[a][YY], fr->v[a][ZZ]);
-            }
-        }
-        else
-        {
-            fprintf(out, "VELOCITYRED\n");
-            for (i = 0; i < nout; i++)
-            {
-                if (index)
-                {
-                    a = index[i];
-                }
-                else
-                {
-                    a = i;
-                }
-                fprintf(out, "%15.9f%15.9f%15.9f\n",
-                        fr->v[a][XX], fr->v[a][YY], fr->v[a][ZZ]);
-            }
-        }
-        fprintf(out, "END\n");
-    }
-    if (fr->bBox)
-    {
-        fprintf(out, "BOX\n");
-        fprintf(out, "%15.9f%15.9f%15.9f",
-                fr->box[XX][XX], fr->box[YY][YY], fr->box[ZZ][ZZ]);
-        if (fr->box[XX][YY] || fr->box[XX][ZZ] || fr->box[YY][XX] ||
-            fr->box[YY][ZZ] || fr->box[ZZ][XX] || fr->box[ZZ][YY])
-        {
-            fprintf(out, "%15.9f%15.9f%15.9f%15.9f%15.9f%15.9f",
-                    fr->box[XX][YY], fr->box[XX][ZZ], fr->box[YY][XX],
-                    fr->box[YY][ZZ], fr->box[ZZ][XX], fr->box[ZZ][YY]);
-        }
-        fprintf(out, "\n");
-        fprintf(out, "END\n");
-    }
-}
-
 static int get_espresso_word(FILE *fp, char word[])
 {
     int  ret, nc, i;
@@ -940,7 +520,7 @@ static gmx_bool get_w_conf(FILE *in, const char *infile, char *title,
     char      *p1, *p2, *p3;
 
     newres  = -1;
-    oldres  = NOTSET; /* Unlikely number for the first residue! */
+    oldres  = -12345; /* Unlikely number for the first residue! */
     ddist   = 0;
 
     /* Read the title and number of atoms */