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, 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, TRUE, 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, TRUE);
 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 // TODO molecule index handling is suspected of being broken here
 241 static void tpx_make_chain_identifiers(t_atoms *atoms, t_block *mols)
 242 {
 243     /* We always assign a new chain number, but save the chain id characters
 244      * for larger molecules.
 245      */
 246     const int chainMinAtoms = 15;
 247
 248     int       chainnum = 0;
 249     char      chainid  = 'A';
 250     bool      outOfIds = false;
 251     for (int m = 0; m < mols->nr; m++)
 252     {
 253         int a0 = mols->index[m];
 254         int a1 = mols->index[m+1];
 255         int c;
 256         if (a1 - a0 >= chainMinAtoms && !outOfIds)
 257         {
 258             /* Set the chain id for the output */
 259             c = chainid;
 260             /* Here we allow for the max possible 2*26+10=62 chain ids */
 261             if (chainid == 'Z')
 262             {
 263                 chainid = 'a';
 264             }
 265             else if (chainid == 'z')
 266             {
 267                 chainid = '0';
 268             }
 269             else if (chainid == '9')
 270             {
 271                 outOfIds = true;
 272             }
 273             else
 274             {
 275                 chainid++;
 276             }
 277         }
 278         else
 279         {
 280             c = ' ';
 281         }
 282         for (int a = a0; a < a1; a++)
 283         {
 284             atoms->resinfo[atoms->atom[a].resind].chainnum = chainnum;
 285             atoms->resinfo[atoms->atom[a].resind].chainid  = c;
 286         }
 287         chainnum++;
 288     }
 289
 290     /* Blank out the chain id if there was only one chain */
 291     if (chainid == 'B')
 292     {
 293         for (int r = 0; r < atoms->nres; r++)
 294         {
 295             atoms->resinfo[r].chainid = ' ';
 296         }
 297     }
 298 }
 299
 300 static void read_stx_conf(const char *infile,
 301                           t_symtab *symtab, char **name, t_atoms *atoms,
 302                           rvec x[], rvec *v, int *ePBC, matrix box)
 303 {
 304     FILE       *in;
 305     t_trxframe  fr;
 306     int         ftp;
 307     char        g96_line[STRLEN+1];
 308
 309     if (atoms->nr == 0)
 310     {
 311         fprintf(stderr, "Warning: Number of atoms in %s is 0\n", infile);
 312     }
 313     else if (atoms->atom == nullptr)
 314     {
 315         gmx_mem("Uninitialized array atom");
 316     }
 317
 318     if (ePBC)
 319     {
 320         *ePBC = -1;
 321     }
 322
 323     ftp = fn2ftp(infile);
 324     switch (ftp)
 325     {
 326         case efGRO:
 327             gmx_gro_read_conf(infile, symtab, name, atoms, x, v, box);
 328             break;
 329         case efG96:
 330             fr.natoms = atoms->nr;
 331             fr.atoms  = atoms;
 332             fr.x      = x;
 333             fr.v      = v;
 334             fr.f      = nullptr;
 335             in        = gmx_fio_fopen(infile, "r");
 336             read_g96_conf(in, infile, name, &fr, symtab, g96_line);
 337             gmx_fio_fclose(in);
 338             copy_mat(fr.box, box);
 339             break;
 340         case efPDB:
 341         case efBRK:
 342         case efENT:
 343             gmx_pdb_read_conf(infile, symtab, name, atoms, x, ePBC, box);
 344             break;
 345         case efESP:
 346             gmx_espresso_read_conf(infile, symtab, name, atoms, x, v, box);
 347             break;
 348         default:
 349             gmx_incons("Not supported in read_stx_conf");
 350     }
 351 }
 352
 353 void readConfAndAtoms(const char *infile,
 354                       t_symtab *symtab, char **name, t_atoms *atoms,
 355                       int *ePBC,
 356                       rvec **x, rvec **v, matrix box)
 357 {
 358     int natoms;
 359     get_stx_coordnum(infile, &natoms);
 360
 361     init_t_atoms(atoms, natoms, (fn2ftp(infile) == efPDB));
 362
 363     bool xIsNull = false;
 364     if (x == nullptr)
 365     {
 366         snew(x, 1);
 367         xIsNull = true;
 368     }
 369     snew(*x, natoms);
 370     if (v)
 371     {
 372         snew(*v, natoms);
 373     }
 374     read_stx_conf(infile,
 375                   symtab, name, atoms,
 376                   *x, (v == nullptr) ? nullptr : *v, ePBC, box);
 377     if (xIsNull)
 378     {
 379         sfree(*x);
 380         sfree(x);
 381     }
 382 }
 383
 384 void readConfAndTopology(const char *infile,
 385                          bool *haveTopology, gmx_mtop_t *mtop,
 386                          int *ePBC,
 387                          rvec **x, rvec **v, matrix box)
 388 {
 389     GMX_RELEASE_ASSERT(mtop != nullptr, "readConfAndTopology requires mtop!=NULL");
 390
 391     if (ePBC != nullptr)
 392     {
 393         *ePBC = -1;
 394     }
 395
 396     *haveTopology = fn2bTPX(infile);
 397     if (*haveTopology)
 398     {
 399         t_tpxheader header;
 400         read_tpxheader(infile, &header, TRUE);
 401         if (x)
 402         {
 403             snew(*x, header.natoms);
 404         }
 405         if (v)
 406         {
 407             snew(*v, header.natoms);
 408         }
 409         int natoms;
 410         int ePBC_tmp
 411             = read_tpx(infile, nullptr, box, &natoms,
 412                        (x == nullptr) ? nullptr : *x, (v == nullptr) ? nullptr : *v, mtop);
 413         if (ePBC != nullptr)
 414         {
 415             *ePBC = ePBC_tmp;
 416         }
 417     }
 418     else
 419     {
 420         t_symtab   symtab;
 421         char      *name;
 422         t_atoms    atoms;
 423
 424         open_symtab(&symtab);
 425
 426         readConfAndAtoms(infile, &symtab, &name, &atoms, ePBC, x, v, box);
 427
 428         init_mtop(mtop);
 429         convertAtomsToMtop(&symtab, put_symtab(&symtab, name), &atoms, mtop);
 430         sfree(name);
 431     }
 432 }
 433
 434 gmx_bool read_tps_conf(const char *infile, t_topology *top, int *ePBC,
 435                        rvec **x, rvec **v, matrix box, gmx_bool requireMasses)
 436 {
 437     bool        haveTopology;
 438     gmx_mtop_t  mtop;
 439
 440     readConfAndTopology(infile, &haveTopology, &mtop, ePBC, x, v, box);
 441
 442     *top = gmx_mtop_t_to_t_topology(&mtop, true);
 443
 444     tpx_make_chain_identifiers(&top->atoms, &top->mols);
 445
 446     if (requireMasses && !top->atoms.haveMass)
 447     {
 448         atomsSetMassesBasedOnNames(&top->atoms, TRUE);
 449
 450         if (!top->atoms.haveMass)
 451         {
 452             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.");
 453         }
 454     }
 455
 456     return haveTopology;
 457 }