src/gromacs/fileio/espio.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2005, The GROMACS development team.
   5  * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 #include "gmxpre.h"
  37
  38 #include "espio.h"
  39
  40 #include <cstdio>
  41 #include <cstdlib>
  42 #include <cstring>
  43
  44 #include "gromacs/fileio/gmxfio.h"
  45 #include "gromacs/math/vec.h"
  46 #include "gromacs/pbcutil/pbc.h"
  47 #include "gromacs/topology/atoms.h"
  48 #include "gromacs/topology/symtab.h"
  49 #include "gromacs/utility/cstringutil.h"
  50 #include "gromacs/utility/fatalerror.h"
  51 #include "gromacs/utility/smalloc.h"
  52
  53 static int get_espresso_word(FILE *fp, char word[])
  54 {
  55     int  ret, nc, i;
  56
  57     ret = 0;
  58     nc  = 0;
  59
  60     do
  61     {
  62         i = fgetc(fp);
  63         if (i != EOF)
  64         {
  65             if (i == ' ' || i == '\n' || i == '\t')
  66             {
  67                 if (nc > 0)
  68                 {
  69                     ret = 1;
  70                 }
  71             }
  72             else if (i == '{')
  73             {
  74                 if (nc == 0)
  75                 {
  76                     word[nc++] = '{';
  77                 }
  78                 ret = 2;
  79             }
  80             else if (i == '}')
  81             {
  82                 if (nc == 0)
  83                 {
  84                     word[nc++] = '}';
  85                 }
  86                 ret = 3;
  87             }
  88             else
  89             {
  90                 word[nc++] = (char)i;
  91             }
  92         }
  93     }
  94     while (i != EOF && ret == 0);
  95
  96     word[nc] = '\0';
  97
  98     return ret;
  99 }
 100
 101 static int check_open_parenthesis(FILE *fp, int r,
 102                                   const char *infile, const char *keyword)
 103 {
 104     int  level_inc;
 105     char word[STRLEN];
 106
 107     level_inc = 0;
 108     if (r == 2)
 109     {
 110         level_inc++;
 111     }
 112     else
 113     {
 114         r = get_espresso_word(fp, word);
 115         if (r == 2)
 116         {
 117             level_inc++;
 118         }
 119         else
 120         {
 121             gmx_fatal(FARGS, "Expected '{' after '%s' in file '%s'",
 122                       keyword, infile);
 123         }
 124     }
 125
 126     return level_inc;
 127 }
 128
 129 static int check_close_parenthesis(FILE *fp, int r,
 130                                    const char *infile, const char *keyword)
 131 {
 132     int  level_inc;
 133     char word[STRLEN];
 134
 135     level_inc = 0;
 136     if (r == 3)
 137     {
 138         level_inc--;
 139     }
 140     else
 141     {
 142         r = get_espresso_word(fp, word);
 143         if (r == 3)
 144         {
 145             level_inc--;
 146         }
 147         else
 148         {
 149             gmx_fatal(FARGS, "Expected '}' after section '%s' in file '%s'",
 150                       keyword, infile);
 151         }
 152     }
 153
 154     return level_inc;
 155 }
 156
 157 enum {
 158     espID, espPOS, espTYPE, espQ, espV, espF, espMOLECULE, espNR
 159 };
 160 static const char *const esp_prop[espNR] = {
 161     "id", "pos", "type", "q", "v", "f",
 162     "molecule"
 163 };
 164
 165 void read_espresso_conf(const char *infile, char *title,
 166                         t_atoms *atoms, rvec x[], rvec *v, matrix box)
 167 {
 168     t_symtab *symtab = NULL;
 169     FILE     *fp;
 170     char      word[STRLEN], buf[STRLEN];
 171     int       level, r, nprop, p, i, m, molnr;
 172     int       prop[32];
 173     double    d;
 174     gmx_bool  bFoundParticles, bFoundProp, bFoundVariable, bMol;
 175
 176     if (!symtab)
 177     {
 178         snew(symtab, 1);
 179         open_symtab(symtab);
 180     }
 181     // TODO: The code does not understand titles it writes...
 182     title[0] = '\0';
 183
 184     clear_mat(box);
 185
 186     fp = gmx_fio_fopen(infile, "r");
 187
 188     bFoundParticles = FALSE;
 189     bFoundVariable  = FALSE;
 190     bMol            = FALSE;
 191     level           = 0;
 192     while ((r = get_espresso_word(fp, word)))
 193     {
 194         if (level == 1 && std::strcmp(word, "particles") == 0 && !bFoundParticles)
 195         {
 196             bFoundParticles = TRUE;
 197             level          += check_open_parenthesis(fp, r, infile, "particles");
 198             nprop           = 0;
 199             while (level == 2 && (r = get_espresso_word(fp, word)))
 200             {
 201                 bFoundProp = FALSE;
 202                 for (p = 0; p < espNR; p++)
 203                 {
 204                     if (strcmp(word, esp_prop[p]) == 0)
 205                     {
 206                         bFoundProp    = TRUE;
 207                         prop[nprop++] = p;
 208                         /* printf("  prop[%d] = %s\n",nprop-1,esp_prop[prop[nprop-1]]); */
 209                     }
 210                 }
 211                 if (!bFoundProp && word[0] != '}')
 212                 {
 213                     gmx_fatal(FARGS, "Can not read Espresso files with particle property '%s'", word);
 214                 }
 215                 if (bFoundProp && p == espMOLECULE)
 216                 {
 217                     bMol = TRUE;
 218                 }
 219                 if (r == 3)
 220                 {
 221                     level--;
 222                 }
 223             }
 224
 225             i = 0;
 226             while (level > 0 && (r = get_espresso_word(fp, word)))
 227             {
 228                 if (r == 2)
 229                 {
 230                     level++;
 231                 }
 232                 else if (r == 3)
 233                 {
 234                     level--;
 235                 }
 236                 if (level == 2)
 237                 {
 238                     for (p = 0; p < nprop; p++)
 239                     {
 240                         switch (prop[p])
 241                         {
 242                             case espID:
 243                                 r = get_espresso_word(fp, word);
 244                                 /* Not used */
 245                                 break;
 246                             case espPOS:
 247                                 for (m = 0; m < 3; m++)
 248                                 {
 249                                     r = get_espresso_word(fp, word);
 250                                     sscanf(word, "%lf", &d);
 251                                     x[i][m] = d;
 252                                 }
 253                                 break;
 254                             case espTYPE:
 255                                 r                   = get_espresso_word(fp, word);
 256                                 atoms->atom[i].type = std::strtol(word, NULL, 10);
 257                                 break;
 258                             case espQ:
 259                                 r = get_espresso_word(fp, word);
 260                                 sscanf(word, "%lf", &d);
 261                                 atoms->atom[i].q = d;
 262                                 break;
 263                             case espV:
 264                                 for (m = 0; m < 3; m++)
 265                                 {
 266                                     r = get_espresso_word(fp, word);
 267                                     sscanf(word, "%lf", &d);
 268                                     v[i][m] = d;
 269                                 }
 270                                 break;
 271                             case espF:
 272                                 for (m = 0; m < 3; m++)
 273                                 {
 274                                     r = get_espresso_word(fp, word);
 275                                     /* not used */
 276                                 }
 277                                 break;
 278                             case espMOLECULE:
 279                                 r     = get_espresso_word(fp, word);
 280                                 molnr = std::strtol(word, NULL, 10);
 281                                 if (i == 0 ||
 282                                     atoms->resinfo[atoms->atom[i-1].resind].nr != molnr)
 283                                 {
 284                                     atoms->atom[i].resind =
 285                                         (i == 0 ? 0 : atoms->atom[i-1].resind+1);
 286                                     atoms->resinfo[atoms->atom[i].resind].nr       = molnr;
 287                                     atoms->resinfo[atoms->atom[i].resind].ic       = ' ';
 288                                     atoms->resinfo[atoms->atom[i].resind].chainid  = ' ';
 289                                     atoms->resinfo[atoms->atom[i].resind].chainnum = molnr; /* Not sure if this is right? */
 290                                 }
 291                                 else
 292                                 {
 293                                     atoms->atom[i].resind = atoms->atom[i-1].resind;
 294                                 }
 295                                 break;
 296                         }
 297                     }
 298                     /* Generate an atom name from the particle type */
 299                     sprintf(buf, "T%d", atoms->atom[i].type);
 300                     atoms->atomname[i] = put_symtab(symtab, buf);
 301                     if (bMol)
 302                     {
 303                         if (i == 0 || atoms->atom[i].resind != atoms->atom[i-1].resind)
 304                         {
 305                             atoms->resinfo[atoms->atom[i].resind].name =
 306                                 put_symtab(symtab, "MOL");
 307                         }
 308                     }
 309                     else
 310                     {
 311                         /* Residue number is the atom number */
 312                         atoms->atom[i].resind = i;
 313                         /* Generate an residue name from the particle type */
 314                         if (atoms->atom[i].type < 26)
 315                         {
 316                             sprintf(buf, "T%c", 'A'+atoms->atom[i].type);
 317                         }
 318                         else
 319                         {
 320                             sprintf(buf, "T%c%c",
 321                                     'A'+atoms->atom[i].type/26, 'A'+atoms->atom[i].type%26);
 322                         }
 323                         t_atoms_set_resinfo(atoms, i, symtab, buf, i, ' ', 0, ' ');
 324                     }
 325
 326                     if (r == 3)
 327                     {
 328                         level--;
 329                     }
 330                     i++;
 331                 }
 332             }
 333             atoms->nres = atoms->nr;
 334
 335             if (i != atoms->nr)
 336             {
 337                 gmx_fatal(FARGS, "Internal inconsistency in Espresso routines, read %d atoms, expected %d atoms", i, atoms->nr);
 338             }
 339         }
 340         else if (level == 1 && std::strcmp(word, "variable") == 0 && !bFoundVariable)
 341         {
 342             bFoundVariable = TRUE;
 343             level         += check_open_parenthesis(fp, r, infile, "variable");
 344             while (level == 2 && (r = get_espresso_word(fp, word)))
 345             {
 346                 if (level == 2 && std::strcmp(word, "box_l") == 0)
 347                 {
 348                     for (m = 0; m < 3; m++)
 349                     {
 350                         r = get_espresso_word(fp, word);
 351                         sscanf(word, "%lf", &d);
 352                         box[m][m] = d;
 353                     }
 354                     level += check_close_parenthesis(fp, r, infile, "box_l");
 355                 }
 356             }
 357         }
 358         else if (r == 2)
 359         {
 360             level++;
 361         }
 362         else if (r == 3)
 363         {
 364             level--;
 365         }
 366     }
 367
 368     if (!bFoundParticles)
 369     {
 370         fprintf(stderr, "Did not find a particles section in Espresso file '%s'\n",
 371                 infile);
 372     }
 373
 374     gmx_fio_fclose(fp);
 375 }
 376
 377 int get_espresso_coordnum(const char *infile)
 378 {
 379     FILE    *fp;
 380     char     word[STRLEN];
 381     int      natoms, level, r;
 382     gmx_bool bFoundParticles;
 383
 384     natoms = 0;
 385
 386     fp = gmx_fio_fopen(infile, "r");
 387
 388     bFoundParticles = FALSE;
 389     level           = 0;
 390     while ((r = get_espresso_word(fp, word)) && !bFoundParticles)
 391     {
 392         if (level == 1 && strcmp(word, "particles") == 0 && !bFoundParticles)
 393         {
 394             bFoundParticles = TRUE;
 395             level          += check_open_parenthesis(fp, r, infile, "particles");
 396             while (level > 0 && (r = get_espresso_word(fp, word)))
 397             {
 398                 if (r == 2)
 399                 {
 400                     level++;
 401                     if (level == 2)
 402                     {
 403                         natoms++;
 404                     }
 405                 }
 406                 else if (r == 3)
 407                 {
 408                     level--;
 409                 }
 410             }
 411         }
 412         else if (r == 2)
 413         {
 414             level++;
 415         }
 416         else if (r == 3)
 417         {
 418             level--;
 419         }
 420     }
 421     if (!bFoundParticles)
 422     {
 423         fprintf(stderr, "Did not find a particles section in Espresso file '%s'\n",
 424                 infile);
 425     }
 426
 427     gmx_fio_fclose(fp);
 428
 429     return natoms;
 430 }
 431
 432 void write_espresso_conf_indexed(FILE *out, const char *title,
 433                                  t_atoms *atoms, int nx, atom_id *index,
 434                                  rvec *x, rvec *v, matrix box)
 435 {
 436     int i, j;
 437
 438     fprintf(out, "# %s\n", title);
 439     if (TRICLINIC(box))
 440     {
 441         gmx_warning("The Espresso format does not support triclinic unit-cells");
 442     }
 443     fprintf(out, "{variable {box_l %f %f %f}}\n", box[0][0], box[1][1], box[2][2]);
 444
 445     fprintf(out, "{particles {id pos type q%s}\n", v ? " v" : "");
 446     for (i = 0; i < nx; i++)
 447     {
 448         if (index)
 449         {
 450             j = index[i];
 451         }
 452         else
 453         {
 454             j = i;
 455         }
 456         fprintf(out, "\t{%d %f %f %f %d %g",
 457                 j, x[j][XX], x[j][YY], x[j][ZZ],
 458                 atoms->atom[j].type, atoms->atom[j].q);
 459         if (v)
 460         {
 461             fprintf(out, " %f %f %f", v[j][XX], v[j][YY], v[j][ZZ]);
 462         }
 463         fprintf(out, "}\n");
 464     }
 465     fprintf(out, "}\n");
 466 }