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
write_sto_conf_indexed(const char * outfile,const char * title,const t_atoms * atoms,const rvec x[],const rvec * v,PbcType pbcType,const matrix box,int nindex,int index[])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 PbcType pbcType,
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, pbcType, box, ' ', -1, nindex, index,
111 nullptr, 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
write_sto_conf(const char * outfile,const char * title,const t_atoms * atoms,const rvec x[],const rvec * v,PbcType pbcType,const matrix box)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 PbcType pbcType,
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, pbcType, 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
write_sto_conf_mtop(const char * outfile,const char * title,const gmx_mtop_t * mtop,const rvec x[],const rvec * v,PbcType pbcType,const matrix box)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 PbcType pbcType,
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, pbcType, box);
202
203 done_atom(&atoms);
204 break;
205 }
206 }
207
get_stx_coordnum(const char * infile,int * natoms)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
fill(t_atoms * atoms,const int startAtomIndex,const int endAtomIndex)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
clearIfNeeded(t_atoms * atoms) const307 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
makeChainIdentifiersAfterTprReading(t_atoms * atoms,const gmx::RangePartitioning & mols)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
tpx_make_chain_identifiers(t_atoms * atoms,const t_block * mols)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
read_stx_conf(const char * infile,t_symtab * symtab,char ** name,t_atoms * atoms,rvec x[],rvec * v,PbcType * pbcType,matrix box)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 PbcType* pbcType,
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 (pbcType)
365 {
366 *pbcType = PbcType::Unset;
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, pbcType, 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
readConfAndAtoms(const char * infile,t_symtab * symtab,char ** name,t_atoms * atoms,PbcType * pbcType,rvec ** x,rvec ** v,matrix box)392 void readConfAndAtoms(const char* infile,
393 t_symtab* symtab,
394 char** name,
395 t_atoms* atoms,
396 PbcType* pbcType,
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, pbcType, 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, pbcType, box);
441 if (xIsNull)
442 {
443 sfree(*x);
444 sfree(x);
445 }
446 }
447
readConfAndTopology(const char * infile,bool * haveTopology,gmx_mtop_t * mtop,PbcType * pbcType,rvec ** x,rvec ** v,matrix box)448 void readConfAndTopology(const char* infile,
449 bool* haveTopology,
450 gmx_mtop_t* mtop,
451 PbcType* pbcType,
452 rvec** x,
453 rvec** v,
454 matrix box)
455 {
456 GMX_RELEASE_ASSERT(mtop != nullptr, "readConfAndTopology requires mtop!=NULL");
457
458 if (pbcType != nullptr)
459 {
460 *pbcType = PbcType::Unset;
461 }
462
463 *haveTopology = fn2bTPX(infile);
464 if (*haveTopology)
465 {
466 TpxFileHeader header = readTpxHeader(infile, true);
467 if (x)
468 {
469 snew(*x, header.natoms);
470 }
471 if (v)
472 {
473 snew(*v, header.natoms);
474 }
475 int natoms;
476 PbcType pbcType_tmp = read_tpx(infile, nullptr, box, &natoms, (x == nullptr) ? nullptr : *x,
477 (v == nullptr) ? nullptr : *v, mtop);
478 if (pbcType != nullptr)
479 {
480 *pbcType = pbcType_tmp;
481 }
482 }
483 else
484 {
485 t_symtab symtab;
486 char* name;
487 t_atoms atoms;
488
489 open_symtab(&symtab);
490
491 readConfAndAtoms(infile, &symtab, &name, &atoms, pbcType, x, v, box);
492
493 convertAtomsToMtop(&symtab, put_symtab(&symtab, name), &atoms, mtop);
494 sfree(name);
495 }
496 }
497
read_tps_conf(const char * infile,t_topology * top,PbcType * pbcType,rvec ** x,rvec ** v,matrix box,gmx_bool requireMasses)498 gmx_bool read_tps_conf(const char* infile, t_topology* top, PbcType* pbcType, rvec** x, rvec** v, matrix box, gmx_bool requireMasses)
499 {
500 bool haveTopology;
501 gmx_mtop_t mtop;
502
503 readConfAndTopology(infile, &haveTopology, &mtop, pbcType, x, v, box);
504
505 *top = gmx_mtop_t_to_t_topology(&mtop, true);
506
507 tpx_make_chain_identifiers(&top->atoms, &top->mols);
508
509 if (requireMasses && !top->atoms.haveMass)
510 {
511 atomsSetMassesBasedOnNames(&top->atoms, TRUE);
512
513 if (!top->atoms.haveMass)
514 {
515 gmx_fatal(FARGS,
516 "Masses were requested, but for some atom(s) masses could not be found in "
517 "the database. Use a tpr file as input, if possible, or add these atoms to "
518 "the mass database.");
519 }
520 }
521
522 return haveTopology;
523 }
524