2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
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"
64 void write_sto_conf_indexed(const char* outfile,
78 ftp = fn2ftp(outfile);
82 out = gmx_fio_fopen(outfile, "w");
83 write_hconf_indexed_p(out, title, atoms, nindex, index, x, v, box);
87 clear_trxframe(&fr, TRUE);
88 fr.natoms = atoms->nr;
90 fr.atoms = const_cast<t_atoms*>(atoms);
92 fr.x = const_cast<rvec*>(x);
96 fr.v = const_cast<rvec*>(v);
99 copy_mat(box, fr.box);
100 out = gmx_fio_fopen(outfile, "w");
101 write_g96_conf(out, title, &fr, nindex, index);
108 out = gmx_fio_fopen(outfile, "w");
109 write_pdbfile_indexed(out, title, atoms, x, ePBC, box, ' ', -1, nindex, index, nullptr,
114 out = gmx_fio_fopen(outfile, "w");
115 write_espresso_conf_indexed(out, title, atoms, nindex, index, x, v, box);
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");
123 void write_sto_conf(const char* outfile,
125 const t_atoms* atoms,
135 ftp = fn2ftp(outfile);
138 case efGRO: write_conf_p(outfile, title, atoms, x, v, box); break;
140 clear_trxframe(&fr, TRUE);
141 fr.natoms = atoms->nr;
143 fr.atoms = const_cast<t_atoms*>(atoms); // TODO check
145 fr.x = const_cast<rvec*>(x);
149 fr.v = const_cast<rvec*>(v);
152 copy_mat(box, fr.box);
153 out = gmx_fio_fopen(outfile, "w");
154 write_g96_conf(out, title, &fr, -1, nullptr);
160 out = gmx_fio_fopen(outfile, "w");
161 write_pdbfile(out, title, atoms, x, ePBC, box, ' ', -1, nullptr);
165 out = gmx_fio_fopen(outfile, "w");
166 write_espresso_conf_indexed(out, title, atoms, atoms->nr, nullptr, x, v, box);
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");
174 void write_sto_conf_mtop(const char* outfile,
176 const gmx_mtop_t* mtop,
186 ftp = fn2ftp(outfile);
190 out = gmx_fio_fopen(outfile, "w");
191 write_hconf_mtop(out, title, mtop, x, v, box);
195 /* This is a brute force approach which requires a lot of memory.
196 * We should implement mtop versions of all writing routines.
198 atoms = gmx_mtop_global_atoms(mtop);
200 write_sto_conf(outfile, title, &atoms, x, v, ePBC, box);
207 static void get_stx_coordnum(const char* infile, int* natoms)
212 char g96_line[STRLEN + 1];
214 ftp = fn2ftp(infile);
215 range_check(ftp, 0, efNR);
218 case efGRO: get_coordnum(infile, natoms); break;
221 in = gmx_fio_fopen(infile, "r");
227 *natoms = read_g96_conf(in, infile, nullptr, &fr, nullptr, g96_line);
234 in = gmx_fio_fopen(infile, "r");
235 get_pdb_coordnum(in, natoms);
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));
243 //! Constructs plausible chain IDs for multi-molecule systems, e.g. when read from .tpr files
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;
253 //! Minimum size for a chain worth giving an ID
254 static constexpr int s_chainMinAtoms = 15;
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;
264 void ChainIdFiller::fill(t_atoms* atoms, const int startAtomIndex, const int endAtomIndex)
266 // TODO remove these some time, extra braces added for review convenience
268 // We always assign a new chain number, but only assign a chain id
269 // characters for larger molecules.
271 if (endAtomIndex - startAtomIndex >= s_chainMinAtoms && !outOfIds_)
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')
280 else if (nextChainId_ == 'z')
284 else if (nextChainId_ == '9')
295 chainIdToAssign = ' ';
297 for (int a = startAtomIndex; a < endAtomIndex; a++)
299 atoms->resinfo[atoms->atom[a].resind].chainnum = nextChainNumber_;
300 atoms->resinfo[atoms->atom[a].resind].chainid = chainIdToAssign;
306 void ChainIdFiller::clearIfNeeded(t_atoms* atoms) const
308 /* Blank out the chain id if there was only one chain */
309 if (nextChainId_ == 'B')
311 for (int r = 0; r < atoms->nres; r++)
313 atoms->resinfo[r].chainid = ' ';
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)
321 ChainIdFiller filler;
322 for (auto m = 0; m != mols.numBlocks(); ++m)
324 filler.fill(atoms, mols.block(m).begin(), mols.block(m).end());
326 filler.clearIfNeeded(atoms);
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)
332 ChainIdFiller filler;
333 for (int m = 0; m < mols->nr; m++)
335 filler.fill(atoms, mols->index[m], mols->index[m + 1]);
337 filler.clearIfNeeded(atoms);
340 static void read_stx_conf(const char* infile,
352 char g96_line[STRLEN + 1];
356 fprintf(stderr, "Warning: Number of atoms in %s is 0\n", infile);
358 else if (atoms->atom == nullptr)
360 gmx_mem("Uninitialized array atom");
368 ftp = fn2ftp(infile);
371 case efGRO: gmx_gro_read_conf(infile, symtab, name, atoms, x, v, box); break;
373 fr.natoms = atoms->nr;
378 in = gmx_fio_fopen(infile, "r");
379 read_g96_conf(in, infile, name, &fr, symtab, g96_line);
381 copy_mat(fr.box, box);
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");
391 void readConfAndAtoms(const char* infile,
400 GMX_RELEASE_ASSERT(infile, "Need a valid file name string");
402 if (fn2ftp(infile) == efTPR)
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);
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++.
417 mtop.symtab.symbuf = nullptr;
424 get_stx_coordnum(infile, &natoms);
426 init_t_atoms(atoms, natoms, (fn2ftp(infile) == efPDB));
428 bool xIsNull = false;
439 read_stx_conf(infile, symtab, name, atoms, *x, (v == nullptr) ? nullptr : *v, ePBC, box);
447 void readConfAndTopology(const char* infile, bool* haveTopology, gmx_mtop_t* mtop, int* ePBC, rvec** x, rvec** v, matrix box)
449 GMX_RELEASE_ASSERT(mtop != nullptr, "readConfAndTopology requires mtop!=NULL");
456 *haveTopology = fn2bTPX(infile);
459 TpxFileHeader header = readTpxHeader(infile, true);
462 snew(*x, header.natoms);
466 snew(*v, header.natoms);
469 int ePBC_tmp = read_tpx(infile, nullptr, box, &natoms, (x == nullptr) ? nullptr : *x,
470 (v == nullptr) ? nullptr : *v, mtop);
482 open_symtab(&symtab);
484 readConfAndAtoms(infile, &symtab, &name, &atoms, ePBC, x, v, box);
486 convertAtomsToMtop(&symtab, put_symtab(&symtab, name), &atoms, mtop);
491 gmx_bool read_tps_conf(const char* infile, t_topology* top, int* ePBC, rvec** x, rvec** v, matrix box, gmx_bool requireMasses)
496 readConfAndTopology(infile, &haveTopology, &mtop, ePBC, x, v, box);
498 *top = gmx_mtop_t_to_t_topology(&mtop, true);
500 tpx_make_chain_identifiers(&top->atoms, &top->mols);
502 if (requireMasses && !top->atoms.haveMass)
504 atomsSetMassesBasedOnNames(&top->atoms, TRUE);
506 if (!top->atoms.haveMass)
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.");