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