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, 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.
47 #include "gromacs/fileio/gmxfio.h"
48 #include "gromacs/topology/atoms.h"
49 #include "gromacs/topology/mtop_util.h"
50 #include "gromacs/topology/symtab.h"
51 #include "gromacs/topology/topology.h"
52 #include "gromacs/trajectory/trajectoryframe.h"
53 #include "gromacs/utility/coolstuff.h"
54 #include "gromacs/utility/cstringutil.h"
55 #include "gromacs/utility/fatalerror.h"
56 #include "gromacs/utility/futil.h"
57 #include "gromacs/utility/smalloc.h"
59 static void get_coordnum_fp(FILE *in
, char *title
, int *natoms
)
63 fgets2(title
, STRLEN
, in
);
64 fgets2(line
, STRLEN
, in
);
65 if (sscanf(line
, "%d", natoms
) != 1)
67 gmx_fatal(FARGS
, "gro file does not have the number of atoms on the second line");
71 void get_coordnum(const char *infile
, int *natoms
)
76 in
= gmx_fio_fopen(infile
, "r");
77 get_coordnum_fp(in
, title
, natoms
);
81 static gmx_bool
get_w_conf(FILE *in
, const char *infile
, char *title
,
82 t_symtab
*symtab
, t_atoms
*atoms
, int *ndec
,
83 rvec x
[], rvec
*v
, matrix box
)
86 char resname
[6], oldresname
[6];
87 char line
[STRLEN
+1], *ptr
;
89 double x1
, y1
, z1
, x2
, y2
, z2
;
91 int natoms
, i
, m
, resnr
, newres
, oldres
, ddist
, c
;
92 gmx_bool bFirst
, bVel
, oldResFirst
;
100 /* Read the title and number of atoms */
101 get_coordnum_fp(in
, title
, &natoms
);
103 if (natoms
> atoms
->nr
)
105 gmx_fatal(FARGS
, "gro file contains more atoms (%d) than expected (%d)",
108 else if (natoms
< atoms
->nr
)
110 fprintf(stderr
, "Warning: gro file contains less atoms (%d) than expected"
111 " (%d)\n", natoms
, atoms
->nr
);
119 oldresname
[0] = '\0';
121 /* just pray the arrays are big enough */
122 for (i
= 0; (i
< natoms
); i
++)
124 if ((fgets2(line
, STRLEN
, in
)) == NULL
)
126 gmx_fatal(FARGS
, "Unexpected end of file in file %s at line %d",
129 if (strlen(line
) < 39)
131 gmx_fatal(FARGS
, "Invalid line in %s for atom %d:\n%s", infile
, i
+1, line
);
134 /* determine read precision from distance between periods
139 p1
= strchr(line
, '.');
142 gmx_fatal(FARGS
, "A coordinate in file %s does not contain a '.'", infile
);
144 p2
= strchr(&p1
[1], '.');
147 gmx_fatal(FARGS
, "A coordinate in file %s does not contain a '.'", infile
);
152 p3
= strchr(&p2
[1], '.');
155 gmx_fatal(FARGS
, "A coordinate in file %s does not contain a '.'", infile
);
158 if (p3
- p2
!= ddist
)
160 gmx_fatal(FARGS
, "The spacing of the decimal points in file %s is not consistent for x, y and z", infile
);
165 memcpy(name
, line
, 5);
167 sscanf(name
, "%d", &resnr
);
168 sscanf(line
+5, "%5s", resname
);
170 if (!oldResFirst
|| oldres
!= resnr
|| strncmp(resname
, oldresname
, sizeof(resname
)))
175 if (newres
>= natoms
)
177 gmx_fatal(FARGS
, "More residues than atoms in %s (natoms = %d)",
180 atoms
->atom
[i
].resind
= newres
;
181 t_atoms_set_resinfo(atoms
, i
, symtab
, resname
, resnr
, ' ', 0, ' ');
185 atoms
->atom
[i
].resind
= newres
;
189 std::memcpy(name
, line
+10, 5);
190 atoms
->atomname
[i
] = put_symtab(symtab
, name
);
192 /* Copy resname to oldresname after we are done with the sanity check above */
193 std::strncpy(oldresname
, resname
, sizeof(oldresname
));
195 /* eventueel controle atomnumber met i+1 */
197 /* coordinates (start after residue data) */
199 /* Read fixed format */
200 for (m
= 0; m
< DIM
; m
++)
202 for (c
= 0; (c
< ddist
&& ptr
[0]); c
++)
208 if (sscanf(buf
, "%lf %lf", &x1
, &x2
) != 1)
210 gmx_fatal(FARGS
, "Something is wrong in the coordinate formatting of file %s. Note that gro is fixed format (see the manual)", infile
);
218 /* velocities (start after residues and coordinates) */
221 /* Read fixed format */
222 for (m
= 0; m
< DIM
; m
++)
224 for (c
= 0; (c
< ddist
&& ptr
[0]); c
++)
230 if (sscanf(buf
, "%lf", &x1
) != 1)
242 atoms
->nres
= newres
+ 1;
245 fgets2(line
, STRLEN
, in
);
246 if (sscanf(line
, "%lf%lf%lf", &x1
, &y1
, &z1
) != 3)
248 gmx_warning("Bad box in file %s", infile
);
250 /* Generate a cubic box */
251 for (m
= 0; (m
< DIM
); m
++)
253 xmin
[m
] = xmax
[m
] = x
[0][m
];
255 for (i
= 1; (i
< atoms
->nr
); i
++)
257 for (m
= 0; (m
< DIM
); m
++)
259 xmin
[m
] = std::min(xmin
[m
], x
[i
][m
]);
260 xmax
[m
] = std::max(xmax
[m
], x
[i
][m
]);
263 for (i
= 0; i
< DIM
; i
++)
265 for (m
= 0; m
< DIM
; m
++)
270 for (m
= 0; (m
< DIM
); m
++)
272 box
[m
][m
] = (xmax
[m
]-xmin
[m
]);
274 fprintf(stderr
, "Generated a cubic box %8.3f x %8.3f x %8.3f\n",
275 box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
]);
279 /* We found the first three values, the diagonal elements */
283 if (sscanf (line
, "%*f%*f%*f%lf%lf%lf%lf%lf%lf",
284 &x1
, &y1
, &z1
, &x2
, &y2
, &z2
) != 6)
286 x1
= y1
= z1
= x2
= y2
= z2
= 0.0;
299 void gmx_gro_read_conf(const char *infile
,
300 t_topology
*top
, rvec x
[], rvec
*v
, matrix box
)
302 FILE *in
= gmx_fio_fopen(infile
, "r");
305 get_w_conf(in
, infile
, title
, &top
->symtab
, &top
->atoms
, &ndec
, x
, v
, box
);
306 top
->name
= put_symtab(&top
->symtab
, title
);
310 static gmx_bool
gmx_one_before_eof(FILE *fp
)
315 if ((beof
= fread(data
, 1, 1, fp
)) == 1)
317 gmx_fseek(fp
, -1, SEEK_CUR
);
322 gmx_bool
gro_next_x_or_v(FILE *status
, t_trxframe
*fr
)
326 char title
[STRLEN
], *p
;
330 if (gmx_one_before_eof(status
))
335 open_symtab(&symtab
);
336 atoms
.nr
= fr
->natoms
;
337 snew(atoms
.atom
, fr
->natoms
);
338 atoms
.nres
= fr
->natoms
;
339 snew(atoms
.resinfo
, fr
->natoms
);
340 snew(atoms
.atomname
, fr
->natoms
);
342 fr
->bV
= get_w_conf(status
, title
, title
, &symtab
, &atoms
, &ndec
, fr
->x
, fr
->v
, fr
->box
);
345 /* prec = 10^ndec: */
346 for (i
= 0; i
< ndec
; i
++)
356 sfree(atoms
.resinfo
);
357 sfree(atoms
.atomname
);
358 done_symtab(&symtab
);
360 if ((p
= strstr(title
, "t=")) != NULL
)
363 if (sscanf(p
, "%lf", &tt
) == 1)
375 if (atoms
.nr
!= fr
->natoms
)
377 gmx_fatal(FARGS
, "Number of atoms in gro frame (%d) doesn't match the number in the previous frame (%d)", atoms
.nr
, fr
->natoms
);
383 int gro_first_x_or_v(FILE *status
, t_trxframe
*fr
)
388 fprintf(stderr
, "Reading frames from gro file");
389 get_coordnum_fp(status
, title
, &fr
->natoms
);
391 fprintf(stderr
, " '%s', %d atoms.\n", title
, fr
->natoms
);
396 gmx_file("No coordinates in gro file");
399 snew(fr
->x
, fr
->natoms
);
400 snew(fr
->v
, fr
->natoms
);
401 gro_next_x_or_v(status
, fr
);
406 static void make_hconf_format(int pr
, gmx_bool bVel
, char format
[])
410 /* build format string for printing,
411 something like "%8.3f" for x and "%8.4f" for v */
424 sprintf(format
, "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
425 l
, pr
, l
, pr
, l
, pr
, l
, vpr
, l
, vpr
, l
, vpr
);
429 sprintf(format
, "%%%d.%df%%%d.%df%%%d.%df\n", l
, pr
, l
, pr
, l
, pr
);
434 static void write_hconf_box(FILE *out
, int pr
, const matrix box
)
445 if (box
[XX
][YY
] || box
[XX
][ZZ
] || box
[YY
][XX
] || box
[YY
][ZZ
] ||
446 box
[ZZ
][XX
] || box
[ZZ
][YY
])
448 sprintf(format
, "%%%d.%df%%%d.%df%%%d.%df"
449 "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
450 l
, pr
, l
, pr
, l
, pr
, l
, pr
, l
, pr
, l
, pr
, l
, pr
, l
, pr
, l
, pr
);
452 box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
],
453 box
[XX
][YY
], box
[XX
][ZZ
], box
[YY
][XX
],
454 box
[YY
][ZZ
], box
[ZZ
][XX
], box
[ZZ
][YY
]);
458 sprintf(format
, "%%%d.%df%%%d.%df%%%d.%df\n", l
, pr
, l
, pr
, l
, pr
);
460 box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
]);
464 void write_hconf_indexed_p(FILE *out
, const char *title
, const t_atoms
*atoms
,
465 int nx
, const int index
[], int pr
,
466 const rvec
*x
, const rvec
*v
, const matrix box
)
468 char resnm
[6], nm
[6], format
[100];
469 int ai
, i
, resind
, resnr
;
471 fprintf(out
, "%s\n", (title
&& title
[0]) ? title
: gmx::bromacs().c_str());
472 fprintf(out
, "%5d\n", nx
);
474 make_hconf_format(pr
, v
!= NULL
, format
);
476 for (i
= 0; (i
< nx
); i
++)
480 resind
= atoms
->atom
[ai
].resind
;
481 std::strncpy(resnm
, " ??? ", sizeof(resnm
)-1);
482 if (resind
< atoms
->nres
)
484 std::strncpy(resnm
, *atoms
->resinfo
[resind
].name
, sizeof(resnm
)-1);
485 resnr
= atoms
->resinfo
[resind
].nr
;
489 std::strncpy(resnm
, " ??? ", sizeof(resnm
)-1);
495 std::strncpy(nm
, *atoms
->atomname
[ai
], sizeof(nm
)-1);
499 std::strncpy(nm
, " ??? ", sizeof(nm
)-1);
502 fprintf(out
, "%5d%-5.5s%5.5s%5d", resnr
%100000, resnm
, nm
, (ai
+1)%100000);
503 /* next fprintf uses built format string */
507 x
[ai
][XX
], x
[ai
][YY
], x
[ai
][ZZ
], v
[ai
][XX
], v
[ai
][YY
], v
[ai
][ZZ
]);
512 x
[ai
][XX
], x
[ai
][YY
], x
[ai
][ZZ
]);
516 write_hconf_box(out
, pr
, box
);
521 void write_hconf_mtop(FILE *out
, const char *title
, gmx_mtop_t
*mtop
, int pr
,
522 const rvec
*x
, const rvec
*v
, const matrix box
)
526 gmx_mtop_atomloop_all_t aloop
;
528 char *atomname
, *resname
;
530 fprintf(out
, "%s\n", (title
&& title
[0]) ? title
: gmx::bromacs().c_str());
531 fprintf(out
, "%5d\n", mtop
->natoms
);
533 make_hconf_format(pr
, v
!= NULL
, format
);
535 aloop
= gmx_mtop_atomloop_all_init(mtop
);
536 while (gmx_mtop_atomloop_all_next(aloop
, &i
, &atom
))
538 gmx_mtop_atomloop_all_names(aloop
, &atomname
, &resnr
, &resname
);
540 fprintf(out
, "%5d%-5.5s%5.5s%5d",
541 resnr
%100000, resname
, atomname
, (i
+1)%100000);
542 /* next fprintf uses built format string */
546 x
[i
][XX
], x
[i
][YY
], x
[i
][ZZ
], v
[i
][XX
], v
[i
][YY
], v
[i
][ZZ
]);
551 x
[i
][XX
], x
[i
][YY
], x
[i
][ZZ
]);
555 write_hconf_box(out
, pr
, box
);
560 void write_hconf_p(FILE *out
, const char *title
, const t_atoms
*atoms
, int pr
,
561 const rvec
*x
, const rvec
*v
, const matrix box
)
567 for (i
= 0; (i
< atoms
->nr
); i
++)
571 write_hconf_indexed_p(out
, title
, atoms
, atoms
->nr
, aa
, pr
, x
, v
, box
);
575 void write_conf_p(const char *outfile
, const char *title
,
576 const t_atoms
*atoms
, int pr
,
577 const rvec
*x
, const rvec
*v
, const matrix box
)
581 out
= gmx_fio_fopen(outfile
, "w");
582 write_hconf_p(out
, title
, atoms
, pr
, x
, v
, box
);