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,2016, 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/topology/topology.h"
  50 #include "gromacs/utility/cstringutil.h"
  51 #include "gromacs/utility/fatalerror.h"
  52 #include "gromacs/utility/smalloc.h"
  53
  54 static int get_espresso_word(FILE *fp, char word[])
  55 {
  56     int  ret, nc, i;
  57
  58     ret = 0;
  59     nc  = 0;
  60
  61     do
  62     {
  63         i = fgetc(fp);
  64         if (i != EOF)
  65         {
  66             if (i == ' ' || i == '\n' || i == '\t')
  67             {
  68                 if (nc > 0)
  69                 {
  70                     ret = 1;
  71                 }
  72             }
  73             else if (i == '{')
  74             {
  75                 if (nc == 0)
  76                 {
  77                     word[nc++] = '{';
  78                 }
  79                 ret = 2;
  80             }
  81             else if (i == '}')
  82             {
  83                 if (nc == 0)
  84                 {
  85                     word[nc++] = '}';
  86                 }
  87                 ret = 3;
  88             }
  89             else
  90             {
  91                 word[nc++] = (char)i;
  92             }
  93         }
  94     }
  95     while (i != EOF && ret == 0);
  96
  97     word[nc] = '\0';
  98
  99     return ret;
 100 }
 101
 102 static int check_open_parenthesis(FILE *fp, int r,
 103                                   const char *infile, const char *keyword)
 104 {
 105     int  level_inc;
 106     char word[STRLEN];
 107
 108     level_inc = 0;
 109     if (r == 2)
 110     {
 111         level_inc++;
 112     }
 113     else
 114     {
 115         r = get_espresso_word(fp, word);
 116         if (r == 2)
 117         {
 118             level_inc++;
 119         }
 120         else
 121         {
 122             gmx_fatal(FARGS, "Expected '{' after '%s' in file '%s'",
 123                       keyword, infile);
 124         }
 125     }
 126
 127     return level_inc;
 128 }
 129
 130 static int check_close_parenthesis(FILE *fp, int r,
 131                                    const char *infile, const char *keyword)
 132 {
 133     int  level_inc;
 134     char word[STRLEN];
 135
 136     level_inc = 0;
 137     if (r == 3)
 138     {
 139         level_inc--;
 140     }
 141     else
 142     {
 143         r = get_espresso_word(fp, word);
 144         if (r == 3)
 145         {
 146             level_inc--;
 147         }
 148         else
 149         {
 150             gmx_fatal(FARGS, "Expected '}' after section '%s' in file '%s'",
 151                       keyword, infile);
 152         }
 153     }
 154
 155     return level_inc;
 156 }
 157
 158 enum {
 159     espID, espPOS, espTYPE, espQ, espV, espF, espMOLECULE, espNR
 160 };
 161 static const char *const esp_prop[espNR] = {
 162     "id", "pos", "type", "q", "v", "f",
 163     "molecule"
 164 };
 165
 166 void gmx_espresso_read_conf(const char *infile,
 167                             t_symtab *symtab, char ***name, t_atoms *atoms,
 168                             rvec x[], rvec *v, matrix box)
 169 {
 170     FILE     *fp;
 171     char      word[STRLEN], buf[STRLEN];
 172     int       level, r, nprop, p, i, m, molnr;
 173     int       prop[32];
 174     double    d;
 175     gmx_bool  bFoundParticles, bFoundProp, bFoundVariable, bMol;
 176
 177     // No title reading implemented for espresso files
 178     *name = put_symtab(symtab, "");
 179
 180     clear_mat(box);
 181
 182     atoms->haveMass    = FALSE;
 183     atoms->haveCharge  = FALSE;
 184     atoms->haveType    = FALSE;
 185     atoms->haveBState  = FALSE;
 186     atoms->havePdbInfo = FALSE;
 187
 188     fp = gmx_fio_fopen(infile, "r");
 189
 190     bFoundParticles = FALSE;
 191     bFoundVariable  = FALSE;
 192     bMol            = FALSE;
 193     level           = 0;
 194     while ((r = get_espresso_word(fp, word)))
 195     {
 196         if (level == 1 && std::strcmp(word, "particles") == 0 && !bFoundParticles)
 197         {
 198             bFoundParticles = TRUE;
 199             level          += check_open_parenthesis(fp, r, infile, "particles");
 200             nprop           = 0;
 201             while (level == 2 && (r = get_espresso_word(fp, word)))
 202             {
 203                 bFoundProp = FALSE;
 204                 for (p = 0; p < espNR; p++)
 205                 {
 206                     if (strcmp(word, esp_prop[p]) == 0)
 207                     {
 208                         bFoundProp    = TRUE;
 209                         prop[nprop++] = p;
 210                         if (p == espQ)
 211                         {
 212                             atoms->haveCharge = TRUE;
 213                         }
 214
 215                         if (debug)
 216                         {
 217                             fprintf(debug, "  prop[%d] = %s\n", nprop-1, esp_prop[prop[nprop-1]]);
 218                         }
 219                     }
 220                 }
 221                 if (!bFoundProp && word[0] != '}')
 222                 {
 223                     gmx_fatal(FARGS, "Can not read Espresso files with particle property '%s'", word);
 224                 }
 225                 if (bFoundProp && p == espMOLECULE)
 226                 {
 227                     bMol = TRUE;
 228                 }
 229                 if (r == 3)
 230                 {
 231                     level--;
 232                 }
 233             }
 234
 235             i = 0;
 236             while (level > 0 && (r = get_espresso_word(fp, word)))
 237             {
 238                 if (r == 2)
 239                 {
 240                     level++;
 241                 }
 242                 else if (r == 3)
 243                 {
 244                     level--;
 245                 }
 246                 if (level == 2)
 247                 {
 248                     for (p = 0; p < nprop; p++)
 249                     {
 250                         switch (prop[p])
 251                         {
 252                             case espID:
 253                                 r = get_espresso_word(fp, word);
 254                                 /* Not used */
 255                                 break;
 256                             case espPOS:
 257                                 for (m = 0; m < 3; m++)
 258                                 {
 259                                     r = get_espresso_word(fp, word);
 260                                     sscanf(word, "%lf", &d);
 261                                     x[i][m] = d;
 262                                 }
 263                                 break;
 264                             case espTYPE:
 265                                 r                   = get_espresso_word(fp, word);
 266                                 atoms->atom[i].type = std::strtol(word, NULL, 10);
 267                                 break;
 268                             case espQ:
 269                                 r = get_espresso_word(fp, word);
 270                                 sscanf(word, "%lf", &d);
 271                                 atoms->atom[i].q = d;
 272                                 break;
 273                             case espV:
 274                                 for (m = 0; m < 3; m++)
 275                                 {
 276                                     r = get_espresso_word(fp, word);
 277                                     sscanf(word, "%lf", &d);
 278                                     v[i][m] = d;
 279                                 }
 280                                 break;
 281                             case espF:
 282                                 for (m = 0; m < 3; m++)
 283                                 {
 284                                     r = get_espresso_word(fp, word);
 285                                     /* not used */
 286                                 }
 287                                 break;
 288                             case espMOLECULE:
 289                                 r     = get_espresso_word(fp, word);
 290                                 molnr = std::strtol(word, NULL, 10);
 291                                 if (i == 0 ||
 292                                     atoms->resinfo[atoms->atom[i-1].resind].nr != molnr)
 293                                 {
 294                                     atoms->atom[i].resind =
 295                                         (i == 0 ? 0 : atoms->atom[i-1].resind+1);
 296                                     atoms->resinfo[atoms->atom[i].resind].nr       = molnr;
 297                                     atoms->resinfo[atoms->atom[i].resind].ic       = ' ';
 298                                     atoms->resinfo[atoms->atom[i].resind].chainid  = ' ';
 299                                     atoms->resinfo[atoms->atom[i].resind].chainnum = molnr; /* Not sure if this is right? */
 300                                 }
 301                                 else
 302                                 {
 303                                     atoms->atom[i].resind = atoms->atom[i-1].resind;
 304                                 }
 305                                 break;
 306                         }
 307                     }
 308                     /* Generate an atom name from the particle type */
 309                     sprintf(buf, "T%hu", atoms->atom[i].type);
 310                     atoms->atomname[i] = put_symtab(symtab, buf);
 311                     if (bMol)
 312                     {
 313                         if (i == 0 || atoms->atom[i].resind != atoms->atom[i-1].resind)
 314                         {
 315                             atoms->resinfo[atoms->atom[i].resind].name =
 316                                 put_symtab(symtab, "MOL");
 317                         }
 318                     }
 319                     else
 320                     {
 321                         /* Residue number is the atom number */
 322                         atoms->atom[i].resind = i;
 323                         /* Generate an residue name from the particle type */
 324                         if (atoms->atom[i].type < 26)
 325                         {
 326                             sprintf(buf, "T%c", 'A'+atoms->atom[i].type);
 327                         }
 328                         else
 329                         {
 330                             sprintf(buf, "T%c%c",
 331                                     'A'+atoms->atom[i].type/26, 'A'+atoms->atom[i].type%26);
 332                         }
 333                         t_atoms_set_resinfo(atoms, i, symtab, buf, i, ' ', 0, ' ');
 334                     }
 335
 336                     if (r == 3)
 337                     {
 338                         level--;
 339                     }
 340                     i++;
 341                 }
 342             }
 343             atoms->nres = atoms->nr;
 344
 345             if (i != atoms->nr)
 346             {
 347                 gmx_fatal(FARGS, "Internal inconsistency in Espresso routines, read %d atoms, expected %d atoms", i, atoms->nr);
 348             }
 349         }
 350         else if (level == 1 && std::strcmp(word, "variable") == 0 && !bFoundVariable)
 351         {
 352             bFoundVariable = TRUE;
 353             level         += check_open_parenthesis(fp, r, infile, "variable");
 354             while (level == 2 && (r = get_espresso_word(fp, word)))
 355             {
 356                 if (level == 2 && std::strcmp(word, "box_l") == 0)
 357                 {
 358                     for (m = 0; m < 3; m++)
 359                     {
 360                         r = get_espresso_word(fp, word);
 361                         sscanf(word, "%lf", &d);
 362                         box[m][m] = d;
 363                     }
 364                     level += check_close_parenthesis(fp, r, infile, "box_l");
 365                 }
 366             }
 367         }
 368         else if (r == 2)
 369         {
 370             level++;
 371         }
 372         else if (r == 3)
 373         {
 374             level--;
 375         }
 376     }
 377
 378     if (!bFoundParticles)
 379     {
 380         fprintf(stderr, "Did not find a particles section in Espresso file '%s'\n",
 381                 infile);
 382     }
 383
 384     gmx_fio_fclose(fp);
 385 }
 386
 387 int get_espresso_coordnum(const char *infile)
 388 {
 389     FILE    *fp;
 390     char     word[STRLEN];
 391     int      natoms, level, r;
 392     gmx_bool bFoundParticles;
 393
 394     natoms = 0;
 395
 396     fp = gmx_fio_fopen(infile, "r");
 397
 398     bFoundParticles = FALSE;
 399     level           = 0;
 400     while ((r = get_espresso_word(fp, word)) && !bFoundParticles)
 401     {
 402         if (level == 1 && strcmp(word, "particles") == 0 && !bFoundParticles)
 403         {
 404             bFoundParticles = TRUE;
 405             level          += check_open_parenthesis(fp, r, infile, "particles");
 406             while (level > 0 && (r = get_espresso_word(fp, word)))
 407             {
 408                 if (r == 2)
 409                 {
 410                     level++;
 411                     if (level == 2)
 412                     {
 413                         natoms++;
 414                     }
 415                 }
 416                 else if (r == 3)
 417                 {
 418                     level--;
 419                 }
 420             }
 421         }
 422         else if (r == 2)
 423         {
 424             level++;
 425         }
 426         else if (r == 3)
 427         {
 428             level--;
 429         }
 430     }
 431     if (!bFoundParticles)
 432     {
 433         fprintf(stderr, "Did not find a particles section in Espresso file '%s'\n",
 434                 infile);
 435     }
 436
 437     gmx_fio_fclose(fp);
 438
 439     return natoms;
 440 }
 441
 442 void write_espresso_conf_indexed(FILE *out, const char *title,
 443                                  const t_atoms *atoms, int nx, int *index,
 444                                  const rvec *x, const rvec *v, const matrix box)
 445 {
 446     int i, j;
 447
 448     fprintf(out, "# %s\n", title);
 449     if (TRICLINIC(box))
 450     {
 451         gmx_warning("The Espresso format does not support triclinic unit-cells");
 452     }
 453     fprintf(out, "{variable {box_l %f %f %f}}\n", box[0][0], box[1][1], box[2][2]);
 454
 455     fprintf(out, "{particles {id pos type q%s}\n", v ? " v" : "");
 456     for (i = 0; i < nx; i++)
 457     {
 458         if (index)
 459         {
 460             j = index[i];
 461         }
 462         else
 463         {
 464             j = i;
 465         }
 466         fprintf(out, "\t{%d %f %f %f %hu %g",
 467                 j, x[j][XX], x[j][YY], x[j][ZZ],
 468                 atoms->atom[j].type, atoms->atom[j].q);
 469         if (v)
 470         {
 471             fprintf(out, " %f %f %f", v[j][XX], v[j][YY], v[j][ZZ]);
 472         }
 473         fprintf(out, "}\n");
 474     }
 475     fprintf(out, "}\n");
 476 }