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, 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.
46 #include "gromacs/fileio/gmxfio.h"
47 #include "gromacs/fileio/trx.h"
48 #include "gromacs/legacyheaders/copyrite.h"
49 #include "gromacs/topology/atoms.h"
50 #include "gromacs/topology/mtop_util.h"
51 #include "gromacs/topology/symtab.h"
52 #include "gromacs/topology/topology.h"
53 #include "gromacs/utility/cstringutil.h"
54 #include "gromacs/utility/fatalerror.h"
55 #include "gromacs/utility/futil.h"
56 #include "gromacs/utility/smalloc.h"
58 static void get_coordnum_fp(FILE *in, char *title, int *natoms)
62 fgets2(title, STRLEN, in);
63 fgets2(line, STRLEN, in);
64 if (sscanf(line, "%d", natoms) != 1)
66 gmx_fatal(FARGS, "gro file does not have the number of atoms on the second line");
70 void get_coordnum(const char *infile, int *natoms)
75 in = gmx_fio_fopen(infile, "r");
76 get_coordnum_fp(in, title, natoms);
80 static gmx_bool get_w_conf(FILE *in, const char *infile, char *title,
81 t_symtab *symtab, t_atoms *atoms, int *ndec,
82 rvec x[], rvec *v, matrix box)
85 char resname[6], oldresname[6];
86 char line[STRLEN+1], *ptr;
88 double x1, y1, z1, x2, y2, z2;
90 int natoms, i, m, resnr, newres, oldres, ddist, c;
91 gmx_bool bFirst, bVel;
95 oldres = -12345; /* Unlikely number for the first residue! */
98 /* Read the title and number of atoms */
99 get_coordnum_fp(in, title, &natoms);
101 if (natoms > atoms->nr)
103 gmx_fatal(FARGS, "gro file contains more atoms (%d) than expected (%d)",
106 else if (natoms < atoms->nr)
108 fprintf(stderr, "Warning: gro file contains less atoms (%d) than expected"
109 " (%d)\n", natoms, atoms->nr);
117 oldresname[0] = '\0';
119 /* just pray the arrays are big enough */
120 for (i = 0; (i < natoms); i++)
122 if ((fgets2(line, STRLEN, in)) == NULL)
124 gmx_fatal(FARGS, "Unexpected end of file in file %s at line %d",
127 if (strlen(line) < 39)
129 gmx_fatal(FARGS, "Invalid line in %s for atom %d:\n%s", infile, i+1, line);
132 /* determine read precision from distance between periods
137 p1 = strchr(line, '.');
140 gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
142 p2 = strchr(&p1[1], '.');
145 gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
150 p3 = strchr(&p2[1], '.');
153 gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
156 if (p3 - p2 != ddist)
158 gmx_fatal(FARGS, "The spacing of the decimal points in file %s is not consistent for x, y and z", infile);
163 memcpy(name, line, 5);
165 sscanf(name, "%d", &resnr);
166 sscanf(line+5, "%5s", resname);
168 if (resnr != oldres || strncmp(resname, oldresname, sizeof(resname)))
172 if (newres >= natoms)
174 gmx_fatal(FARGS, "More residues than atoms in %s (natoms = %d)",
177 atoms->atom[i].resind = newres;
178 t_atoms_set_resinfo(atoms, i, symtab, resname, resnr, ' ', 0, ' ');
182 atoms->atom[i].resind = newres;
186 std::memcpy(name, line+10, 5);
187 atoms->atomname[i] = put_symtab(symtab, name);
189 /* Copy resname to oldresname after we are done with the sanity check above */
190 std::strncpy(oldresname, resname, sizeof(oldresname));
192 /* eventueel controle atomnumber met i+1 */
194 /* coordinates (start after residue data) */
196 /* Read fixed format */
197 for (m = 0; m < DIM; m++)
199 for (c = 0; (c < ddist && ptr[0]); c++)
205 if (sscanf(buf, "%lf %lf", &x1, &x2) != 1)
207 gmx_fatal(FARGS, "Something is wrong in the coordinate formatting of file %s. Note that gro is fixed format (see the manual)", infile);
215 /* velocities (start after residues and coordinates) */
218 /* Read fixed format */
219 for (m = 0; m < DIM; m++)
221 for (c = 0; (c < ddist && ptr[0]); c++)
227 if (sscanf(buf, "%lf", &x1) != 1)
239 atoms->nres = newres + 1;
242 fgets2(line, STRLEN, in);
243 if (sscanf(line, "%lf%lf%lf", &x1, &y1, &z1) != 3)
245 gmx_warning("Bad box in file %s", infile);
247 /* Generate a cubic box */
248 for (m = 0; (m < DIM); m++)
250 xmin[m] = xmax[m] = x[0][m];
252 for (i = 1; (i < atoms->nr); i++)
254 for (m = 0; (m < DIM); m++)
256 xmin[m] = std::min(xmin[m], x[i][m]);
257 xmax[m] = std::max(xmax[m], x[i][m]);
260 for (i = 0; i < DIM; i++)
262 for (m = 0; m < DIM; m++)
267 for (m = 0; (m < DIM); m++)
269 box[m][m] = (xmax[m]-xmin[m]);
271 fprintf(stderr, "Generated a cubic box %8.3f x %8.3f x %8.3f\n",
272 box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
276 /* We found the first three values, the diagonal elements */
280 if (sscanf (line, "%*f%*f%*f%lf%lf%lf%lf%lf%lf",
281 &x1, &y1, &z1, &x2, &y2, &z2) != 6)
283 x1 = y1 = z1 = x2 = y2 = z2 = 0.0;
296 void read_whole_conf(const char *infile, char *title,
297 t_atoms *atoms, rvec x[], rvec *v, matrix box)
304 in = gmx_fio_fopen(infile, "r");
306 open_symtab(&symtab);
307 get_w_conf(in, infile, title, &symtab, atoms, &ndec, x, v, box);
308 /* We can't free the symbols, as they are still used in atoms, so
309 * the only choice is to leak them. */
310 free_symtab(&symtab);
315 static gmx_bool gmx_one_before_eof(FILE *fp)
320 if ((beof = fread(data, 1, 1, fp)) == 1)
322 gmx_fseek(fp, -1, SEEK_CUR);
327 gmx_bool gro_next_x_or_v(FILE *status, t_trxframe *fr)
331 char title[STRLEN], *p;
335 if (gmx_one_before_eof(status))
340 open_symtab(&symtab);
341 atoms.nr = fr->natoms;
342 snew(atoms.atom, fr->natoms);
343 atoms.nres = fr->natoms;
344 snew(atoms.resinfo, fr->natoms);
345 snew(atoms.atomname, fr->natoms);
347 fr->bV = get_w_conf(status, title, title, &symtab, &atoms, &ndec, fr->x, fr->v, fr->box);
350 /* prec = 10^ndec: */
351 for (i = 0; i < ndec; i++)
361 sfree(atoms.resinfo);
362 sfree(atoms.atomname);
363 done_symtab(&symtab);
365 if ((p = strstr(title, "t=")) != NULL)
368 if (sscanf(p, "%lf", &tt) == 1)
380 if (atoms.nr != fr->natoms)
382 gmx_fatal(FARGS, "Number of atoms in gro frame (%d) doesn't match the number in the previous frame (%d)", atoms.nr, fr->natoms);
388 int gro_first_x_or_v(FILE *status, t_trxframe *fr)
393 fprintf(stderr, "Reading frames from gro file");
394 get_coordnum_fp(status, title, &fr->natoms);
396 fprintf(stderr, " '%s', %d atoms.\n", title, fr->natoms);
401 gmx_file("No coordinates in gro file");
404 snew(fr->x, fr->natoms);
405 snew(fr->v, fr->natoms);
406 gro_next_x_or_v(status, fr);
411 static void make_hconf_format(int pr, gmx_bool bVel, char format[])
415 /* build format string for printing,
416 something like "%8.3f" for x and "%8.4f" for v */
429 sprintf(format, "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
430 l, pr, l, pr, l, pr, l, vpr, l, vpr, l, vpr);
434 sprintf(format, "%%%d.%df%%%d.%df%%%d.%df\n", l, pr, l, pr, l, pr);
439 static void write_hconf_box(FILE *out, int pr, matrix box)
450 if (box[XX][YY] || box[XX][ZZ] || box[YY][XX] || box[YY][ZZ] ||
451 box[ZZ][XX] || box[ZZ][YY])
453 sprintf(format, "%%%d.%df%%%d.%df%%%d.%df"
454 "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
455 l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr);
457 box[XX][XX], box[YY][YY], box[ZZ][ZZ],
458 box[XX][YY], box[XX][ZZ], box[YY][XX],
459 box[YY][ZZ], box[ZZ][XX], box[ZZ][YY]);
463 sprintf(format, "%%%d.%df%%%d.%df%%%d.%df\n", l, pr, l, pr, l, pr);
465 box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
469 void write_hconf_indexed_p(FILE *out, const char *title, t_atoms *atoms,
470 int nx, const atom_id index[], int pr,
471 rvec *x, rvec *v, matrix box)
473 char resnm[6], nm[6], format[100];
474 int ai, i, resind, resnr;
477 fprintf(out, "%s\n", (title && title[0]) ? title : format);
478 fprintf(out, "%5d\n", nx);
480 make_hconf_format(pr, v != NULL, format);
482 for (i = 0; (i < nx); i++)
486 resind = atoms->atom[ai].resind;
487 std::strncpy(resnm, " ??? ", sizeof(resnm)-1);
488 if (resind < atoms->nres)
490 std::strncpy(resnm, *atoms->resinfo[resind].name, sizeof(resnm)-1);
491 resnr = atoms->resinfo[resind].nr;
495 std::strncpy(resnm, " ??? ", sizeof(resnm)-1);
501 std::strncpy(nm, *atoms->atomname[ai], sizeof(nm)-1);
505 std::strncpy(nm, " ??? ", sizeof(nm)-1);
508 fprintf(out, "%5d%-5.5s%5.5s%5d", resnr%100000, resnm, nm, (ai+1)%100000);
509 /* next fprintf uses built format string */
513 x[ai][XX], x[ai][YY], x[ai][ZZ], v[ai][XX], v[ai][YY], v[ai][ZZ]);
518 x[ai][XX], x[ai][YY], x[ai][ZZ]);
522 write_hconf_box(out, pr, box);
527 void write_hconf_mtop(FILE *out, const char *title, gmx_mtop_t *mtop, int pr,
528 rvec *x, rvec *v, matrix box)
532 gmx_mtop_atomloop_all_t aloop;
534 char *atomname, *resname;
537 fprintf(out, "%s\n", (title && title[0]) ? title : format);
538 fprintf(out, "%5d\n", mtop->natoms);
540 make_hconf_format(pr, v != NULL, format);
542 aloop = gmx_mtop_atomloop_all_init(mtop);
543 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
545 gmx_mtop_atomloop_all_names(aloop, &atomname, &resnr, &resname);
547 fprintf(out, "%5d%-5.5s%5.5s%5d",
548 resnr%100000, resname, atomname, (i+1)%100000);
549 /* next fprintf uses built format string */
553 x[i][XX], x[i][YY], x[i][ZZ], v[i][XX], v[i][YY], v[i][ZZ]);
558 x[i][XX], x[i][YY], x[i][ZZ]);
562 write_hconf_box(out, pr, box);
567 void write_hconf_p(FILE *out, const char *title, t_atoms *atoms, int pr,
568 rvec *x, rvec *v, matrix box)
574 for (i = 0; (i < atoms->nr); i++)
578 write_hconf_indexed_p(out, title, atoms, atoms->nr, aa, pr, x, v, box);
582 void write_conf_p(const char *outfile, const char *title,
583 t_atoms *atoms, int pr,
584 rvec *x, rvec *v, matrix box)
588 out = gmx_fio_fopen(outfile, "w");
589 write_hconf_p(out, title, atoms, pr, x, v, box);