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