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,
  65                             const char*    title,
  66                             const t_atoms* atoms,
  67                             const rvec     x[],
  68                             const rvec*    v,
  69                             int            ePBC,
  70                             const matrix   box,
  71                             int            nindex,
  72                             int            index[])
  73 {
  74     FILE*      out;
  75     int        ftp;
  76     t_trxframe fr;
  77
  78     ftp = fn2ftp(outfile);
  79     switch (ftp)
  80     {
  81         case efGRO:
  82             out = gmx_fio_fopen(outfile, "w");
  83             write_hconf_indexed_p(out, title, atoms, nindex, index, x, v, box);
  84             gmx_fio_fclose(out);
  85             break;
  86         case efG96:
  87             clear_trxframe(&fr, TRUE);
  88             fr.natoms = atoms->nr;
  89             fr.bAtoms = TRUE;
  90             fr.atoms  = const_cast<t_atoms*>(atoms);
  91             fr.bX     = TRUE;
  92             fr.x      = const_cast<rvec*>(x);
  93             if (v)
  94             {
  95                 fr.bV = TRUE;
  96                 fr.v  = const_cast<rvec*>(v);
  97             }
  98             fr.bBox = TRUE;
  99             copy_mat(box, fr.box);
 100             out = gmx_fio_fopen(outfile, "w");
 101             write_g96_conf(out, title, &fr, nindex, index);
 102             gmx_fio_fclose(out);
 103             break;
 104         case efPDB:
 105         case efBRK:
 106         case efENT:
 107         case efPQR:
 108             out = gmx_fio_fopen(outfile, "w");
 109             write_pdbfile_indexed(out, title, atoms, x, ePBC, box, ' ', -1, nindex, index, nullptr,
 110                                   ftp == efPQR);
 111             gmx_fio_fclose(out);
 112             break;
 113         case efESP:
 114             out = gmx_fio_fopen(outfile, "w");
 115             write_espresso_conf_indexed(out, title, atoms, nindex, index, x, v, box);
 116             gmx_fio_fclose(out);
 117             break;
 118         case efTPR: gmx_fatal(FARGS, "Sorry, can not write a topology to %s", outfile);
 119         default: gmx_incons("Not supported in write_sto_conf_indexed");
 120     }
 121 }
 122
 123 void write_sto_conf(const char*    outfile,
 124                     const char*    title,
 125                     const t_atoms* atoms,
 126                     const rvec     x[],
 127                     const rvec*    v,
 128                     int            ePBC,
 129                     const matrix   box)
 130 {
 131     FILE*      out;
 132     int        ftp;
 133     t_trxframe fr;
 134
 135     ftp = fn2ftp(outfile);
 136     switch (ftp)
 137     {
 138         case efGRO: write_conf_p(outfile, title, atoms, x, v, box); break;
 139         case efG96:
 140             clear_trxframe(&fr, TRUE);
 141             fr.natoms = atoms->nr;
 142             fr.bAtoms = TRUE;
 143             fr.atoms  = const_cast<t_atoms*>(atoms); // TODO check
 144             fr.bX     = TRUE;
 145             fr.x      = const_cast<rvec*>(x);
 146             if (v)
 147             {
 148                 fr.bV = TRUE;
 149                 fr.v  = const_cast<rvec*>(v);
 150             }
 151             fr.bBox = TRUE;
 152             copy_mat(box, fr.box);
 153             out = gmx_fio_fopen(outfile, "w");
 154             write_g96_conf(out, title, &fr, -1, nullptr);
 155             gmx_fio_fclose(out);
 156             break;
 157         case efPDB:
 158         case efBRK:
 159         case efENT:
 160             out = gmx_fio_fopen(outfile, "w");
 161             write_pdbfile(out, title, atoms, x, ePBC, box, ' ', -1, nullptr);
 162             gmx_fio_fclose(out);
 163             break;
 164         case efESP:
 165             out = gmx_fio_fopen(outfile, "w");
 166             write_espresso_conf_indexed(out, title, atoms, atoms->nr, nullptr, x, v, box);
 167             gmx_fio_fclose(out);
 168             break;
 169         case efTPR: gmx_fatal(FARGS, "Sorry, can not write a topology to %s", outfile);
 170         default: gmx_incons("Not supported in write_sto_conf");
 171     }
 172 }
 173
 174 void write_sto_conf_mtop(const char*       outfile,
 175                          const char*       title,
 176                          const gmx_mtop_t* mtop,
 177                          const rvec        x[],
 178                          const rvec*       v,
 179                          int               ePBC,
 180                          const matrix      box)
 181 {
 182     int     ftp;
 183     FILE*   out;
 184     t_atoms atoms;
 185
 186     ftp = fn2ftp(outfile);
 187     switch (ftp)
 188     {
 189         case efGRO:
 190             out = gmx_fio_fopen(outfile, "w");
 191             write_hconf_mtop(out, title, mtop, x, v, box);
 192             gmx_fio_fclose(out);
 193             break;
 194         default:
 195             /* This is a brute force approach which requires a lot of memory.
 196              * We should implement mtop versions of all writing routines.
 197              */
 198             atoms = gmx_mtop_global_atoms(mtop);
 199
 200             write_sto_conf(outfile, title, &atoms, x, v, ePBC, box);
 201
 202             done_atom(&atoms);
 203             break;
 204     }
 205 }
 206
 207 static void get_stx_coordnum(const char* infile, int* natoms)
 208 {
 209     FILE*      in;
 210     int        ftp;
 211     t_trxframe fr;
 212     char       g96_line[STRLEN + 1];
 213
 214     ftp = fn2ftp(infile);
 215     range_check(ftp, 0, efNR);
 216     switch (ftp)
 217     {
 218         case efGRO: get_coordnum(infile, natoms); break;
 219         case efG96:
 220         {
 221             in        = gmx_fio_fopen(infile, "r");
 222             fr.natoms = -1;
 223             fr.atoms  = nullptr;
 224             fr.x      = nullptr;
 225             fr.v      = nullptr;
 226             fr.f      = nullptr;
 227             *natoms   = read_g96_conf(in, infile, nullptr, &fr, nullptr, g96_line);
 228             gmx_fio_fclose(in);
 229             break;
 230         }
 231         case efPDB:
 232         case efBRK:
 233         case efENT:
 234             in = gmx_fio_fopen(infile, "r");
 235             get_pdb_coordnum(in, natoms);
 236             gmx_fio_fclose(in);
 237             break;
 238         case efESP: *natoms = get_espresso_coordnum(infile); break;
 239         default: gmx_fatal(FARGS, "File type %s not supported in get_stx_coordnum", ftp2ext(ftp));
 240     }
 241 }
 242
 243 //! Constructs plausible chain IDs for multi-molecule systems, e.g. when read from .tpr files
 244 class ChainIdFiller
 245 {
 246 public:
 247     //! Fill in the chain ID for the indicated atom range, which might be a molecule.
 248     void fill(t_atoms* atoms, int startAtomIndex, int endAtomIndex);
 249     //! If only one chain was found, we don't add a chain ID.
 250     void clearIfNeeded(t_atoms* atoms) const;
 251
 252 private:
 253     //! Minimum size for a chain worth giving an ID
 254     static constexpr int s_chainMinAtoms = 15;
 255
 256     //! The number of the next chain that will be assigned.
 257     int nextChainNumber_ = 0;
 258     //! The chain ID of the next chain that will be assigned.
 259     char nextChainId_ = 'A';
 260     //! Whether the set of chain IDs (ie. upper- and lower-case letters and single digits) is exhausted.
 261     bool outOfIds_ = false;
 262 };
 263
 264 void ChainIdFiller::fill(t_atoms* atoms, const int startAtomIndex, const int endAtomIndex)
 265 {
 266     // TODO remove these some time, extra braces added for review convenience
 267     {
 268         // We always assign a new chain number, but only assign a chain id
 269         // characters for larger molecules.
 270         int chainIdToAssign;
 271         if (endAtomIndex - startAtomIndex >= s_chainMinAtoms && !outOfIds_)
 272         {
 273             /* Set the chain id for the output */
 274             chainIdToAssign = nextChainId_;
 275             /* Here we allow for the max possible 2*26+10=62 chain ids */
 276             if (nextChainId_ == 'Z')
 277             {
 278                 nextChainId_ = 'a';
 279             }
 280             else if (nextChainId_ == 'z')
 281             {
 282                 nextChainId_ = '0';
 283             }
 284             else if (nextChainId_ == '9')
 285             {
 286                 outOfIds_ = true;
 287             }
 288             else
 289             {
 290                 nextChainId_++;
 291             }
 292         }
 293         else
 294         {
 295             chainIdToAssign = ' ';
 296         }
 297         for (int a = startAtomIndex; a < endAtomIndex; a++)
 298         {
 299             atoms->resinfo[atoms->atom[a].resind].chainnum = nextChainNumber_;
 300             atoms->resinfo[atoms->atom[a].resind].chainid  = chainIdToAssign;
 301         }
 302         nextChainNumber_++;
 303     }
 304 }
 305
 306 void ChainIdFiller::clearIfNeeded(t_atoms* atoms) const
 307 {
 308     /* Blank out the chain id if there was only one chain */
 309     if (nextChainId_ == 'B')
 310     {
 311         for (int r = 0; r < atoms->nres; r++)
 312         {
 313             atoms->resinfo[r].chainid = ' ';
 314         }
 315     }
 316 }
 317
 318 //! Make chain IDs in the t_atoms for a gmx_mtop_t built from a .tpr file
 319 static void makeChainIdentifiersAfterTprReading(t_atoms* atoms, const gmx::RangePartitioning& mols)
 320 {
 321     ChainIdFiller filler;
 322     for (auto m = 0; m != mols.numBlocks(); ++m)
 323     {
 324         filler.fill(atoms, mols.block(m).begin(), mols.block(m).end());
 325     }
 326     filler.clearIfNeeded(atoms);
 327 }
 328
 329 //! Make chain IDs in the t_atoms for a legacy t_topology built from a .tpr file
 330 static void tpx_make_chain_identifiers(t_atoms* atoms, const t_block* mols)
 331 {
 332     ChainIdFiller filler;
 333     for (int m = 0; m < mols->nr; m++)
 334     {
 335         filler.fill(atoms, mols->index[m], mols->index[m + 1]);
 336     }
 337     filler.clearIfNeeded(atoms);
 338 }
 339
 340 static void read_stx_conf(const char* infile,
 341                           t_symtab*   symtab,
 342                           char**      name,
 343                           t_atoms*    atoms,
 344                           rvec        x[],
 345                           rvec*       v,
 346                           int*        ePBC,
 347                           matrix      box)
 348 {
 349     FILE*      in;
 350     t_trxframe fr;
 351     int        ftp;
 352     char       g96_line[STRLEN + 1];
 353
 354     if (atoms->nr == 0)
 355     {
 356         fprintf(stderr, "Warning: Number of atoms in %s is 0\n", infile);
 357     }
 358     else if (atoms->atom == nullptr)
 359     {
 360         gmx_mem("Uninitialized array atom");
 361     }
 362
 363     if (ePBC)
 364     {
 365         *ePBC = -1;
 366     }
 367
 368     ftp = fn2ftp(infile);
 369     switch (ftp)
 370     {
 371         case efGRO: gmx_gro_read_conf(infile, symtab, name, atoms, x, v, box); 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: gmx_pdb_read_conf(infile, symtab, name, atoms, x, ePBC, box); break;
 386         case efESP: gmx_espresso_read_conf(infile, symtab, name, atoms, x, v, box); break;
 387         default: gmx_incons("Not supported in read_stx_conf");
 388     }
 389 }
 390
 391 void readConfAndAtoms(const char* infile,
 392                       t_symtab*   symtab,
 393                       char**      name,
 394                       t_atoms*    atoms,
 395                       int*        ePBC,
 396                       rvec**      x,
 397                       rvec**      v,
 398                       matrix      box)
 399 {
 400     GMX_RELEASE_ASSERT(infile, "Need a valid file name string");
 401
 402     if (fn2ftp(infile) == efTPR)
 403     {
 404         bool       haveTopology;
 405         gmx_mtop_t mtop;
 406         readConfAndTopology(infile, &haveTopology, &mtop, ePBC, x, v, box);
 407         *symtab                          = mtop.symtab;
 408         *name                            = gmx_strdup(*mtop.name);
 409         *atoms                           = gmx_mtop_global_atoms(&mtop);
 410         gmx::RangePartitioning molecules = gmx_mtop_molecules(mtop);
 411         makeChainIdentifiersAfterTprReading(atoms, molecules);
 412
 413         /* Inelegant solution to avoid all char pointers in atoms becoming
 414          * invalid after destruction of mtop.
 415          * This will be fixed soon by converting t_symtab to C++.
 416          */
 417         mtop.symtab.symbuf = nullptr;
 418         mtop.symtab.nr     = 0;
 419
 420         return;
 421     }
 422
 423     int natoms;
 424     get_stx_coordnum(infile, &natoms);
 425
 426     init_t_atoms(atoms, natoms, (fn2ftp(infile) == efPDB));
 427
 428     bool xIsNull = false;
 429     if (x == nullptr)
 430     {
 431         snew(x, 1);
 432         xIsNull = true;
 433     }
 434     snew(*x, natoms);
 435     if (v)
 436     {
 437         snew(*v, natoms);
 438     }
 439     read_stx_conf(infile, symtab, name, atoms, *x, (v == nullptr) ? nullptr : *v, ePBC, box);
 440     if (xIsNull)
 441     {
 442         sfree(*x);
 443         sfree(x);
 444     }
 445 }
 446
 447 void readConfAndTopology(const char* infile, bool* haveTopology, gmx_mtop_t* mtop, int* ePBC, rvec** x, rvec** v, matrix box)
 448 {
 449     GMX_RELEASE_ASSERT(mtop != nullptr, "readConfAndTopology requires mtop!=NULL");
 450
 451     if (ePBC != nullptr)
 452     {
 453         *ePBC = -1;
 454     }
 455
 456     *haveTopology = fn2bTPX(infile);
 457     if (*haveTopology)
 458     {
 459         TpxFileHeader header = readTpxHeader(infile, true);
 460         if (x)
 461         {
 462             snew(*x, header.natoms);
 463         }
 464         if (v)
 465         {
 466             snew(*v, header.natoms);
 467         }
 468         int natoms;
 469         int ePBC_tmp = read_tpx(infile, nullptr, box, &natoms, (x == nullptr) ? nullptr : *x,
 470                                 (v == nullptr) ? nullptr : *v, mtop);
 471         if (ePBC != nullptr)
 472         {
 473             *ePBC = ePBC_tmp;
 474         }
 475     }
 476     else
 477     {
 478         t_symtab symtab;
 479         char*    name;
 480         t_atoms  atoms;
 481
 482         open_symtab(&symtab);
 483
 484         readConfAndAtoms(infile, &symtab, &name, &atoms, ePBC, x, v, box);
 485
 486         convertAtomsToMtop(&symtab, put_symtab(&symtab, name), &atoms, mtop);
 487         sfree(name);
 488     }
 489 }
 490
 491 gmx_bool read_tps_conf(const char* infile, t_topology* top, int* ePBC, rvec** x, rvec** v, matrix box, gmx_bool requireMasses)
 492 {
 493     bool       haveTopology;
 494     gmx_mtop_t mtop;
 495
 496     readConfAndTopology(infile, &haveTopology, &mtop, ePBC, x, v, box);
 497
 498     *top = gmx_mtop_t_to_t_topology(&mtop, true);
 499
 500     tpx_make_chain_identifiers(&top->atoms, &top->mols);
 501
 502     if (requireMasses && !top->atoms.haveMass)
 503     {
 504         atomsSetMassesBasedOnNames(&top->atoms, TRUE);
 505
 506         if (!top->atoms.haveMass)
 507         {
 508             gmx_fatal(FARGS,
 509                       "Masses were requested, but for some atom(s) masses could not be found in "
 510                       "the database. Use a tpr file as input, if possible, or add these atoms to "
 511                       "the mass database.");
 512         }
 513     }
 514
 515     return haveTopology;
 516 }