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, TRUE);
 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             break;
 115         default:
 116             gmx_incons("Not supported in write_sto_conf_indexed");
 117     }
 118 }
 119
 120 void write_sto_conf(const char *outfile, const char *title, const t_atoms *atoms,
 121                     const rvec x[], const rvec *v, int ePBC, const matrix box)
 122 {
 123     FILE       *out;
 124     int         ftp;
 125     t_trxframe  fr;
 126
 127     ftp = fn2ftp(outfile);
 128     switch (ftp)
 129     {
 130         case efGRO:
 131             write_conf_p(outfile, title, atoms, x, v, box);
 132             break;
 133         case efG96:
 134             clear_trxframe(&fr, TRUE);
 135             fr.natoms = atoms->nr;
 136             fr.bAtoms = TRUE;
 137             fr.atoms  = const_cast<t_atoms *>(atoms); // TODO check
 138             fr.bX     = TRUE;
 139             fr.x      = const_cast<rvec *>(x);
 140             if (v)
 141             {
 142                 fr.bV = TRUE;
 143                 fr.v  = const_cast<rvec *>(v);
 144             }
 145             fr.bBox = TRUE;
 146             copy_mat(box, fr.box);
 147             out = gmx_fio_fopen(outfile, "w");
 148             write_g96_conf(out, title, &fr, -1, nullptr);
 149             gmx_fio_fclose(out);
 150             break;
 151         case efPDB:
 152         case efBRK:
 153         case efENT:
 154             out = gmx_fio_fopen(outfile, "w");
 155             write_pdbfile(out, title, atoms, x, ePBC, box, ' ', -1, nullptr, TRUE);
 156             gmx_fio_fclose(out);
 157             break;
 158         case efESP:
 159             out = gmx_fio_fopen(outfile, "w");
 160             write_espresso_conf_indexed(out, title, atoms, atoms->nr, nullptr, x, v, box);
 161             gmx_fio_fclose(out);
 162             break;
 163         case efTPR:
 164             gmx_fatal(FARGS, "Sorry, can not write a topology to %s", outfile);
 165             break;
 166         default:
 167             gmx_incons("Not supported in write_sto_conf");
 168     }
 169 }
 170
 171 void write_sto_conf_mtop(const char *outfile, const char *title,
 172                          gmx_mtop_t *mtop,
 173                          const rvec x[], const rvec *v, int ePBC, const matrix box)
 174 {
 175     int     ftp;
 176     FILE   *out;
 177     t_atoms atoms;
 178
 179     ftp = fn2ftp(outfile);
 180     switch (ftp)
 181     {
 182         case efGRO:
 183             out = gmx_fio_fopen(outfile, "w");
 184             write_hconf_mtop(out, title, mtop, x, v, box);
 185             gmx_fio_fclose(out);
 186             break;
 187         default:
 188             /* This is a brute force approach which requires a lot of memory.
 189              * We should implement mtop versions of all writing routines.
 190              */
 191             atoms = gmx_mtop_global_atoms(mtop);
 192
 193             write_sto_conf(outfile, title, &atoms, x, v, ePBC, box);
 194
 195             done_atom(&atoms);
 196             break;
 197     }
 198 }
 199
 200 static void get_stx_coordnum(const char *infile, int *natoms)
 201 {
 202     FILE      *in;
 203     int        ftp;
 204     t_trxframe fr;
 205     char       g96_line[STRLEN+1];
 206
 207     ftp = fn2ftp(infile);
 208     range_check(ftp, 0, efNR);
 209     switch (ftp)
 210     {
 211         case efGRO:
 212             get_coordnum(infile, natoms);
 213             break;
 214         case efG96:
 215         {
 216             in        = gmx_fio_fopen(infile, "r");
 217             fr.natoms = -1;
 218             fr.atoms  = nullptr;
 219             fr.x      = nullptr;
 220             fr.v      = nullptr;
 221             fr.f      = nullptr;
 222             *natoms   = read_g96_conf(in, infile, nullptr, &fr, nullptr, g96_line);
 223             gmx_fio_fclose(in);
 224             break;
 225         }
 226         case efPDB:
 227         case efBRK:
 228         case efENT:
 229             in = gmx_fio_fopen(infile, "r");
 230             get_pdb_coordnum(in, natoms);
 231             gmx_fio_fclose(in);
 232             break;
 233         case efESP:
 234             *natoms = get_espresso_coordnum(infile);
 235             break;
 236         default:
 237             gmx_fatal(FARGS, "File type %s not supported in get_stx_coordnum",
 238                       ftp2ext(ftp));
 239     }
 240 }
 241
 242 static void tpx_make_chain_identifiers(t_atoms *atoms, t_block *mols)
 243 {
 244     /* We always assign a new chain number, but save the chain id characters
 245      * for larger molecules.
 246      */
 247     const int chainMinAtoms = 15;
 248
 249     int       chainnum = 0;
 250     char      chainid  = 'A';
 251     bool      outOfIds = false;
 252     for (int m = 0; m < mols->nr; m++)
 253     {
 254         int a0 = mols->index[m];
 255         int a1 = mols->index[m+1];
 256         int c;
 257         if (a1 - a0 >= chainMinAtoms && !outOfIds)
 258         {
 259             /* Set the chain id for the output */
 260             c = chainid;
 261             /* Here we allow for the max possible 2*26+10=62 chain ids */
 262             if (chainid == 'Z')
 263             {
 264                 chainid = 'a';
 265             }
 266             else if (chainid == 'z')
 267             {
 268                 chainid = '0';
 269             }
 270             else if (chainid == '9')
 271             {
 272                 outOfIds = true;
 273             }
 274             else
 275             {
 276                 chainid++;
 277             }
 278         }
 279         else
 280         {
 281             c = ' ';
 282         }
 283         for (int a = a0; a < a1; a++)
 284         {
 285             atoms->resinfo[atoms->atom[a].resind].chainnum = chainnum;
 286             atoms->resinfo[atoms->atom[a].resind].chainid  = c;
 287         }
 288         chainnum++;
 289     }
 290
 291     /* Blank out the chain id if there was only one chain */
 292     if (chainid == 'B')
 293     {
 294         for (int r = 0; r < atoms->nres; r++)
 295         {
 296             atoms->resinfo[r].chainid = ' ';
 297         }
 298     }
 299 }
 300
 301 static void read_stx_conf(const char *infile,
 302                           t_symtab *symtab, char **name, t_atoms *atoms,
 303                           rvec x[], rvec *v, int *ePBC, matrix box)
 304 {
 305     FILE       *in;
 306     t_trxframe  fr;
 307     int         ftp;
 308     char        g96_line[STRLEN+1];
 309
 310     if (atoms->nr == 0)
 311     {
 312         fprintf(stderr, "Warning: Number of atoms in %s is 0\n", infile);
 313     }
 314     else if (atoms->atom == nullptr)
 315     {
 316         gmx_mem("Uninitialized array atom");
 317     }
 318
 319     if (ePBC)
 320     {
 321         *ePBC = -1;
 322     }
 323
 324     ftp = fn2ftp(infile);
 325     switch (ftp)
 326     {
 327         case efGRO:
 328             gmx_gro_read_conf(infile, symtab, name, atoms, x, v, box);
 329             break;
 330         case efG96:
 331             fr.natoms = atoms->nr;
 332             fr.atoms  = atoms;
 333             fr.x      = x;
 334             fr.v      = v;
 335             fr.f      = nullptr;
 336             in        = gmx_fio_fopen(infile, "r");
 337             read_g96_conf(in, infile, name, &fr, symtab, g96_line);
 338             gmx_fio_fclose(in);
 339             copy_mat(fr.box, box);
 340             break;
 341         case efPDB:
 342         case efBRK:
 343         case efENT:
 344             gmx_pdb_read_conf(infile, symtab, name, atoms, x, ePBC, box);
 345             break;
 346         case efESP:
 347             gmx_espresso_read_conf(infile, symtab, name, atoms, x, v, box);
 348             break;
 349         default:
 350             gmx_incons("Not supported in read_stx_conf");
 351     }
 352 }
 353
 354 static void readConfAndAtoms(const char *infile,
 355                              t_symtab *symtab, char **name, t_atoms *atoms,
 356                              int *ePBC,
 357                              rvec **x, rvec **v, matrix box)
 358 {
 359     int natoms;
 360     get_stx_coordnum(infile, &natoms);
 361
 362     init_t_atoms(atoms, natoms, (fn2ftp(infile) == efPDB));
 363
 364     bool xIsNull = false;
 365     if (x == nullptr)
 366     {
 367         snew(x, 1);
 368         xIsNull = true;
 369     }
 370     snew(*x, natoms);
 371     if (v)
 372     {
 373         snew(*v, natoms);
 374     }
 375     read_stx_conf(infile,
 376                   symtab, name, atoms,
 377                   *x, (v == nullptr) ? nullptr : *v, ePBC, box);
 378     if (xIsNull)
 379     {
 380         sfree(*x);
 381         sfree(x);
 382     }
 383 }
 384
 385 void readConfAndTopology(const char *infile,
 386                          bool *haveTopology, gmx_mtop_t *mtop,
 387                          int *ePBC,
 388                          rvec **x, rvec **v, matrix box)
 389 {
 390     GMX_RELEASE_ASSERT(mtop != nullptr, "readConfAndTopology requires mtop!=NULL");
 391
 392     if (ePBC != nullptr)
 393     {
 394         *ePBC = -1;
 395     }
 396
 397     *haveTopology = fn2bTPX(infile);
 398     if (*haveTopology)
 399     {
 400         t_tpxheader header;
 401         read_tpxheader(infile, &header, TRUE);
 402         if (x)
 403         {
 404             snew(*x, header.natoms);
 405         }
 406         if (v)
 407         {
 408             snew(*v, header.natoms);
 409         }
 410         int natoms;
 411         int ePBC_tmp
 412             = read_tpx(infile, nullptr, box, &natoms,
 413                        (x == nullptr) ? nullptr : *x, (v == nullptr) ? nullptr : *v, mtop);
 414         if (ePBC != nullptr)
 415         {
 416             *ePBC = ePBC_tmp;
 417         }
 418     }
 419     else
 420     {
 421         t_symtab   symtab;
 422         char      *name;
 423         t_atoms    atoms;
 424
 425         open_symtab(&symtab);
 426
 427         readConfAndAtoms(infile, &symtab, &name, &atoms, ePBC, x, v, box);
 428
 429         init_mtop(mtop);
 430         convertAtomsToMtop(&symtab, put_symtab(&symtab, name), &atoms, mtop);
 431         sfree(name);
 432     }
 433 }
 434
 435 gmx_bool read_tps_conf(const char *infile, t_topology *top, int *ePBC,
 436                        rvec **x, rvec **v, matrix box, gmx_bool requireMasses)
 437 {
 438     bool        haveTopology;
 439     gmx_mtop_t *mtop;
 440
 441     // Note: We should have an initializer instead of relying on snew
 442     snew(mtop, 1);
 443     readConfAndTopology(infile, &haveTopology, mtop, ePBC, x, v, box);
 444
 445     *top = gmx_mtop_t_to_t_topology(mtop, true);
 446     sfree(mtop);
 447
 448     tpx_make_chain_identifiers(&top->atoms, &top->mols);
 449
 450     if (requireMasses && !top->atoms.haveMass)
 451     {
 452         atomsSetMassesBasedOnNames(&top->atoms, TRUE);
 453
 454         if (!top->atoms.haveMass)
 455         {
 456             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.");
 457         }
 458     }
 459
 460     return haveTopology;
 461 }