src/gromacs/fileio/confio.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "confio.h"
  40
  41 #include <cstdio>
  42 #include <cstring>
  43
  44 #include "gromacs/fileio/espio.h"
  45 #include "gromacs/fileio/filetypes.h"
  46 #include "gromacs/fileio/g96io.h"
  47 #include "gromacs/fileio/gmxfio.h"
  48 #include "gromacs/fileio/groio.h"
  49 #include "gromacs/fileio/pdbio.h"
  50 #include "gromacs/fileio/tpxio.h"
  51 #include "gromacs/fileio/trxio.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/topology/atoms.h"
  54 #include "gromacs/topology/block.h"
  55 #include "gromacs/topology/mtop_util.h"
  56 #include "gromacs/topology/symtab.h"
  57 #include "gromacs/topology/topology.h"
  58 #include "gromacs/trajectory/trajectoryframe.h"
  59 #include "gromacs/utility/cstringutil.h"
  60 #include "gromacs/utility/fatalerror.h"
  61 #include "gromacs/utility/gmxassert.h"
  62 #include "gromacs/utility/smalloc.h"
  63
  64 void write_sto_conf_indexed(const char *outfile, const char *title,
  65                             const t_atoms *atoms,
  66                             const rvec x[], const rvec *v, int ePBC, const matrix box,
  67                             int nindex, int index[])
  68 {
  69     FILE       *out;
  70     int         ftp;
  71     t_trxframe  fr;
  72
  73     ftp = fn2ftp(outfile);
  74     switch (ftp)
  75     {
  76         case efGRO:
  77             out = gmx_fio_fopen(outfile, "w");
  78             write_hconf_indexed_p(out, title, atoms, nindex, index, x, v, box);
  79             gmx_fio_fclose(out);
  80             break;
  81         case efG96:
  82             clear_trxframe(&fr, TRUE);
  83             fr.natoms = atoms->nr;
  84             fr.bAtoms = TRUE;
  85             fr.atoms  = const_cast<t_atoms *>(atoms);
  86             fr.bX     = TRUE;
  87             fr.x      = const_cast<rvec *>(x);
  88             if (v)
  89             {
  90                 fr.bV = TRUE;
  91                 fr.v  = const_cast<rvec *>(v);
  92             }
  93             fr.bBox = TRUE;
  94             copy_mat(box, fr.box);
  95             out = gmx_fio_fopen(outfile, "w");
  96             write_g96_conf(out, title, &fr, nindex, index);
  97             gmx_fio_fclose(out);
  98             break;
  99         case efPDB:
 100         case efBRK:
 101         case efENT:
 102         case efPQR:
 103             out = gmx_fio_fopen(outfile, "w");
 104             write_pdbfile_indexed(out, title, atoms, x, ePBC, box, ' ', -1, nindex, index, nullptr, ftp == efPQR);
 105             gmx_fio_fclose(out);
 106             break;
 107         case efESP:
 108             out = gmx_fio_fopen(outfile, "w");
 109             write_espresso_conf_indexed(out, title, atoms, nindex, index, x, v, box);
 110             gmx_fio_fclose(out);
 111             break;
 112         case efTPR:
 113             gmx_fatal(FARGS, "Sorry, can not write a topology to %s", outfile);
 114         default:
 115             gmx_incons("Not supported in write_sto_conf_indexed");
 116     }
 117 }
 118
 119 void write_sto_conf(const char *outfile, const char *title, const t_atoms *atoms,
 120                     const rvec x[], const rvec *v, int ePBC, const matrix box)
 121 {
 122     FILE       *out;
 123     int         ftp;
 124     t_trxframe  fr;
 125
 126     ftp = fn2ftp(outfile);
 127     switch (ftp)
 128     {
 129         case efGRO:
 130             write_conf_p(outfile, title, atoms, x, v, box);
 131             break;
 132         case efG96:
 133             clear_trxframe(&fr, TRUE);
 134             fr.natoms = atoms->nr;
 135             fr.bAtoms = TRUE;
 136             fr.atoms  = const_cast<t_atoms *>(atoms); // TODO check
 137             fr.bX     = TRUE;
 138             fr.x      = const_cast<rvec *>(x);
 139             if (v)
 140             {
 141                 fr.bV = TRUE;
 142                 fr.v  = const_cast<rvec *>(v);
 143             }
 144             fr.bBox = TRUE;
 145             copy_mat(box, fr.box);
 146             out = gmx_fio_fopen(outfile, "w");
 147             write_g96_conf(out, title, &fr, -1, nullptr);
 148             gmx_fio_fclose(out);
 149             break;
 150         case efPDB:
 151         case efBRK:
 152         case efENT:
 153             out = gmx_fio_fopen(outfile, "w");
 154             write_pdbfile(out, title, atoms, x, ePBC, box, ' ', -1, nullptr);
 155             gmx_fio_fclose(out);
 156             break;
 157         case efESP:
 158             out = gmx_fio_fopen(outfile, "w");
 159             write_espresso_conf_indexed(out, title, atoms, atoms->nr, nullptr, x, v, box);
 160             gmx_fio_fclose(out);
 161             break;
 162         case efTPR:
 163             gmx_fatal(FARGS, "Sorry, can not write a topology to %s", outfile);
 164         default:
 165             gmx_incons("Not supported in write_sto_conf");
 166     }
 167 }
 168
 169 void write_sto_conf_mtop(const char *outfile, const char *title,
 170                          gmx_mtop_t *mtop,
 171                          const rvec x[], const rvec *v, int ePBC, const matrix box)
 172 {
 173     int     ftp;
 174     FILE   *out;
 175     t_atoms atoms;
 176
 177     ftp = fn2ftp(outfile);
 178     switch (ftp)
 179     {
 180         case efGRO:
 181             out = gmx_fio_fopen(outfile, "w");
 182             write_hconf_mtop(out, title, mtop, x, v, box);
 183             gmx_fio_fclose(out);
 184             break;
 185         default:
 186             /* This is a brute force approach which requires a lot of memory.
 187              * We should implement mtop versions of all writing routines.
 188              */
 189             atoms = gmx_mtop_global_atoms(mtop);
 190
 191             write_sto_conf(outfile, title, &atoms, x, v, ePBC, box);
 192
 193             done_atom(&atoms);
 194             break;
 195     }
 196 }
 197
 198 static void get_stx_coordnum(const char *infile, int *natoms)
 199 {
 200     FILE      *in;
 201     int        ftp;
 202     t_trxframe fr;
 203     char       g96_line[STRLEN+1];
 204
 205     ftp = fn2ftp(infile);
 206     range_check(ftp, 0, efNR);
 207     switch (ftp)
 208     {
 209         case efGRO:
 210             get_coordnum(infile, natoms);
 211             break;
 212         case efG96:
 213         {
 214             in        = gmx_fio_fopen(infile, "r");
 215             fr.natoms = -1;
 216             fr.atoms  = nullptr;
 217             fr.x      = nullptr;
 218             fr.v      = nullptr;
 219             fr.f      = nullptr;
 220             *natoms   = read_g96_conf(in, infile, nullptr, &fr, nullptr, g96_line);
 221             gmx_fio_fclose(in);
 222             break;
 223         }
 224         case efPDB:
 225         case efBRK:
 226         case efENT:
 227             in = gmx_fio_fopen(infile, "r");
 228             get_pdb_coordnum(in, natoms);
 229             gmx_fio_fclose(in);
 230             break;
 231         case efESP:
 232             *natoms = get_espresso_coordnum(infile);
 233             break;
 234         default:
 235             gmx_fatal(FARGS, "File type %s not supported in get_stx_coordnum",
 236                       ftp2ext(ftp));
 237     }
 238 }
 239
 240 //! Constructs plausible chain IDs for multi-molecule systems, e.g. when read from .tpr files
 241 class ChainIdFiller
 242 {
 243     public:
 244         //! Fill in the chain ID for the indicated atom range, which might be a molecule.
 245         void fill(t_atoms  *atoms,
 246                   int       startAtomIndex,
 247                   int       endAtomIndex);
 248         //! If only one chain was found, we don't add a chain ID.
 249         void clearIfNeeded(t_atoms *atoms) const;
 250
 251     private:
 252         //! Minimum size for a chain worth giving an ID
 253         static constexpr int s_chainMinAtoms = 15;
 254
 255         //! The number of the next chain that will be assigned.
 256         int  nextChainNumber_ = 0;
 257         //! The chain ID of the next chain that will be assigned.
 258         char nextChainId_ = 'A';
 259         //! Whether the set of chain IDs (ie. upper- and lower-case letters and single digits) is exhausted.
 260         bool outOfIds_ = false;
 261 };
 262
 263 void
 264 ChainIdFiller::fill(t_atoms  *atoms,
 265                     const int startAtomIndex,
 266                     const int endAtomIndex)
 267 {
 268     // TODO remove these some time, extra braces added for review convenience
 269     {
 270         // We always assign a new chain number, but only assign a chain id
 271         // characters for larger molecules.
 272         int chainIdToAssign;
 273         if (endAtomIndex - startAtomIndex >= s_chainMinAtoms && !outOfIds_)
 274         {
 275             /* Set the chain id for the output */
 276             chainIdToAssign = nextChainId_;
 277             /* Here we allow for the max possible 2*26+10=62 chain ids */
 278             if (nextChainId_ == 'Z')
 279             {
 280                 nextChainId_ = 'a';
 281             }
 282             else if (nextChainId_ == 'z')
 283             {
 284                 nextChainId_ = '0';
 285             }
 286             else if (nextChainId_ == '9')
 287             {
 288                 outOfIds_ = true;
 289             }
 290             else
 291             {
 292                 nextChainId_++;
 293             }
 294         }
 295         else
 296         {
 297             chainIdToAssign = ' ';
 298         }
 299         for (int a = startAtomIndex; a < endAtomIndex; a++)
 300         {
 301             atoms->resinfo[atoms->atom[a].resind].chainnum = nextChainNumber_;
 302             atoms->resinfo[atoms->atom[a].resind].chainid  = chainIdToAssign;
 303         }
 304         nextChainNumber_++;
 305     }
 306 }
 307
 308 void
 309 ChainIdFiller::clearIfNeeded(t_atoms *atoms) const
 310 {
 311     /* Blank out the chain id if there was only one chain */
 312     if (nextChainId_ == 'B')
 313     {
 314         for (int r = 0; r < atoms->nres; r++)
 315         {
 316             atoms->resinfo[r].chainid = ' ';
 317         }
 318     }
 319 }
 320
 321 //! Make chain IDs in the t_atoms for a gmx_mtop_t built from a .tpr file
 322 static void makeChainIdentifiersAfterTprReading(t_atoms *atoms, const gmx::RangePartitioning &mols)
 323 {
 324     ChainIdFiller filler;
 325     for (auto m = 0; m != mols.numBlocks(); ++m)
 326     {
 327         filler.fill(atoms, mols.block(m).begin(), mols.block(m).end());
 328     }
 329     filler.clearIfNeeded(atoms);
 330 }
 331
 332 //! Make chain IDs in the t_atoms for a legacy t_topology built from a .tpr file
 333 static void tpx_make_chain_identifiers(t_atoms *atoms, const t_block *mols)
 334 {
 335     ChainIdFiller filler;
 336     for (int m = 0; m < mols->nr; m++)
 337     {
 338         filler.fill(atoms, mols->index[m], mols->index[m+1]);
 339     }
 340     filler.clearIfNeeded(atoms);
 341 }
 342
 343 static void read_stx_conf(const char *infile,
 344                           t_symtab *symtab, char **name, t_atoms *atoms,
 345                           rvec x[], rvec *v, int *ePBC, matrix box)
 346 {
 347     FILE       *in;
 348     t_trxframe  fr;
 349     int         ftp;
 350     char        g96_line[STRLEN+1];
 351
 352     if (atoms->nr == 0)
 353     {
 354         fprintf(stderr, "Warning: Number of atoms in %s is 0\n", infile);
 355     }
 356     else if (atoms->atom == nullptr)
 357     {
 358         gmx_mem("Uninitialized array atom");
 359     }
 360
 361     if (ePBC)
 362     {
 363         *ePBC = -1;
 364     }
 365
 366     ftp = fn2ftp(infile);
 367     switch (ftp)
 368     {
 369         case efGRO:
 370             gmx_gro_read_conf(infile, symtab, name, atoms, x, v, box);
 371             break;
 372         case efG96:
 373             fr.natoms = atoms->nr;
 374             fr.atoms  = atoms;
 375             fr.x      = x;
 376             fr.v      = v;
 377             fr.f      = nullptr;
 378             in        = gmx_fio_fopen(infile, "r");
 379             read_g96_conf(in, infile, name, &fr, symtab, g96_line);
 380             gmx_fio_fclose(in);
 381             copy_mat(fr.box, box);
 382             break;
 383         case efPDB:
 384         case efBRK:
 385         case efENT:
 386             gmx_pdb_read_conf(infile, symtab, name, atoms, x, ePBC, box);
 387             break;
 388         case efESP:
 389             gmx_espresso_read_conf(infile, symtab, name, atoms, x, v, box);
 390             break;
 391         default:
 392             gmx_incons("Not supported in read_stx_conf");
 393     }
 394 }
 395
 396 void readConfAndAtoms(const char *infile,
 397                       t_symtab *symtab, char **name, t_atoms *atoms,
 398                       int *ePBC,
 399                       rvec **x, rvec **v, matrix box)
 400 {
 401     GMX_RELEASE_ASSERT(infile, "Need a valid file name string");
 402
 403     if (fn2ftp(infile) == efTPR)
 404     {
 405         bool       haveTopology;
 406         gmx_mtop_t mtop;
 407         readConfAndTopology(infile, &haveTopology, &mtop, ePBC, x, v, box);
 408         *symtab = mtop.symtab;
 409         *name   = gmx_strdup(*mtop.name);
 410         *atoms  = gmx_mtop_global_atoms(&mtop);
 411         gmx::RangePartitioning molecules = gmx_mtop_molecules(mtop);
 412         makeChainIdentifiersAfterTprReading(atoms, molecules);
 413
 414         /* Inelegant solution to avoid all char pointers in atoms becoming
 415          * invalid after destruction of mtop.
 416          * This will be fixed soon by converting t_symtab to C++.
 417          */
 418         mtop.symtab.symbuf = nullptr;
 419         mtop.symtab.nr     = 0;
 420
 421         return;
 422     }
 423
 424     int natoms;
 425     get_stx_coordnum(infile, &natoms);
 426
 427     init_t_atoms(atoms, natoms, (fn2ftp(infile) == efPDB));
 428
 429     bool xIsNull = false;
 430     if (x == nullptr)
 431     {
 432         snew(x, 1);
 433         xIsNull = true;
 434     }
 435     snew(*x, natoms);
 436     if (v)
 437     {
 438         snew(*v, natoms);
 439     }
 440     read_stx_conf(infile,
 441                   symtab, name, atoms,
 442                   *x, (v == nullptr) ? nullptr : *v, ePBC, box);
 443     if (xIsNull)
 444     {
 445         sfree(*x);
 446         sfree(x);
 447     }
 448 }
 449
 450 void readConfAndTopology(const char *infile,
 451                          bool *haveTopology, gmx_mtop_t *mtop,
 452                          int *ePBC,
 453                          rvec **x, rvec **v, matrix box)
 454 {
 455     GMX_RELEASE_ASSERT(mtop != nullptr, "readConfAndTopology requires mtop!=NULL");
 456
 457     if (ePBC != nullptr)
 458     {
 459         *ePBC = -1;
 460     }
 461
 462     *haveTopology = fn2bTPX(infile);
 463     if (*haveTopology)
 464     {
 465         t_tpxheader header;
 466         read_tpxheader(infile, &header, TRUE);
 467         if (x)
 468         {
 469             snew(*x, header.natoms);
 470         }
 471         if (v)
 472         {
 473             snew(*v, header.natoms);
 474         }
 475         int natoms;
 476         int ePBC_tmp
 477             = read_tpx(infile, nullptr, box, &natoms,
 478                        (x == nullptr) ? nullptr : *x, (v == nullptr) ? nullptr : *v, mtop);
 479         if (ePBC != nullptr)
 480         {
 481             *ePBC = ePBC_tmp;
 482         }
 483     }
 484     else
 485     {
 486         t_symtab   symtab;
 487         char      *name;
 488         t_atoms    atoms;
 489
 490         open_symtab(&symtab);
 491
 492         readConfAndAtoms(infile, &symtab, &name, &atoms, ePBC, x, v, box);
 493
 494         convertAtomsToMtop(&symtab, put_symtab(&symtab, name), &atoms, mtop);
 495         sfree(name);
 496     }
 497 }
 498
 499 gmx_bool read_tps_conf(const char *infile, t_topology *top, int *ePBC,
 500                        rvec **x, rvec **v, matrix box, gmx_bool requireMasses)
 501 {
 502     bool        haveTopology;
 503     gmx_mtop_t  mtop;
 504
 505     readConfAndTopology(infile, &haveTopology, &mtop, ePBC, x, v, box);
 506
 507     *top = gmx_mtop_t_to_t_topology(&mtop, true);
 508
 509     tpx_make_chain_identifiers(&top->atoms, &top->mols);
 510
 511     if (requireMasses && !top->atoms.haveMass)
 512     {
 513         atomsSetMassesBasedOnNames(&top->atoms, TRUE);
 514
 515         if (!top->atoms.haveMass)
 516         {
 517             gmx_fatal(FARGS, "Masses were requested, but for some atom(s) masses could not be found in the database. Use a tpr file as input, if possible, or add these atoms to the mass database.");
 518         }
 519     }
 520
 521     return haveTopology;
 522 }