Split espresso I/O routines from confio.cpp
[alexxy/gromacs.git] / src / gromacs / fileio / confio.cpp
index 14093492279b5a84534198799be21ea7e873b7e3..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"
@@ -55,7 +51,6 @@
 #include "gromacs/fileio/trx.h"
 #include "gromacs/fileio/trxio.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");
-}
-
 void write_sto_conf_indexed(const char *outfile, const char *title,
                             t_atoms *atoms,
                             rvec x[], rvec *v, int ePBC, matrix box,