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