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 /* Note that the .gro reading routine still support variable precision
82 * for backward compatibility with old .gro files.
83 * We have removed writing of variable precision to avoid compatibility
84 * issues with other software packages.
86 static gmx_bool
get_w_conf(FILE *in
, const char *infile
, char *title
,
87 t_symtab
*symtab
, t_atoms
*atoms
, int *ndec
,
88 rvec x
[], rvec
*v
, matrix box
)
91 char resname
[6], oldresname
[6];
92 char line
[STRLEN
+1], *ptr
;
94 double x1
, y1
, z1
, x2
, y2
, z2
;
96 int natoms
, i
, m
, resnr
, newres
, oldres
, ddist
, c
;
97 gmx_bool bFirst
, bVel
, oldResFirst
;
105 /* Read the title and number of atoms */
106 get_coordnum_fp(in
, title
, &natoms
);
108 if (natoms
> atoms
->nr
)
110 gmx_fatal(FARGS
, "gro file contains more atoms (%d) than expected (%d)",
113 else if (natoms
< atoms
->nr
)
115 fprintf(stderr
, "Warning: gro file contains less atoms (%d) than expected"
116 " (%d)\n", natoms
, atoms
->nr
);
119 atoms
->haveMass
= FALSE
;
120 atoms
->haveCharge
= FALSE
;
121 atoms
->haveType
= FALSE
;
122 atoms
->haveBState
= FALSE
;
123 atoms
->havePdbInfo
= FALSE
;
130 oldresname
[0] = '\0';
132 /* just pray the arrays are big enough */
133 for (i
= 0; (i
< natoms
); i
++)
135 if ((fgets2(line
, STRLEN
, in
)) == NULL
)
137 gmx_fatal(FARGS
, "Unexpected end of file in file %s at line %d",
140 if (strlen(line
) < 39)
142 gmx_fatal(FARGS
, "Invalid line in %s for atom %d:\n%s", infile
, i
+1, line
);
145 /* determine read precision from distance between periods
150 p1
= strchr(line
, '.');
153 gmx_fatal(FARGS
, "A coordinate in file %s does not contain a '.'", infile
);
155 p2
= strchr(&p1
[1], '.');
158 gmx_fatal(FARGS
, "A coordinate in file %s does not contain a '.'", infile
);
163 p3
= strchr(&p2
[1], '.');
166 gmx_fatal(FARGS
, "A coordinate in file %s does not contain a '.'", infile
);
169 if (p3
- p2
!= ddist
)
171 gmx_fatal(FARGS
, "The spacing of the decimal points in file %s is not consistent for x, y and z", infile
);
176 memcpy(name
, line
, 5);
178 sscanf(name
, "%d", &resnr
);
179 sscanf(line
+5, "%5s", resname
);
181 if (!oldResFirst
|| oldres
!= resnr
|| strncmp(resname
, oldresname
, sizeof(resname
)))
186 if (newres
>= natoms
)
188 gmx_fatal(FARGS
, "More residues than atoms in %s (natoms = %d)",
191 atoms
->atom
[i
].resind
= newres
;
192 t_atoms_set_resinfo(atoms
, i
, symtab
, resname
, resnr
, ' ', 0, ' ');
196 atoms
->atom
[i
].resind
= newres
;
200 std::memcpy(name
, line
+10, 5);
201 atoms
->atomname
[i
] = put_symtab(symtab
, name
);
203 /* Copy resname to oldresname after we are done with the sanity check above */
204 std::strncpy(oldresname
, resname
, sizeof(oldresname
));
206 /* eventueel controle atomnumber met i+1 */
208 /* coordinates (start after residue data) */
210 /* Read fixed format */
211 for (m
= 0; m
< DIM
; m
++)
213 for (c
= 0; (c
< ddist
&& ptr
[0]); c
++)
219 if (sscanf(buf
, "%lf %lf", &x1
, &x2
) != 1)
221 gmx_fatal(FARGS
, "Something is wrong in the coordinate formatting of file %s. Note that gro is fixed format (see the manual)", infile
);
229 /* velocities (start after residues and coordinates) */
232 /* Read fixed format */
233 for (m
= 0; m
< DIM
; m
++)
235 for (c
= 0; (c
< ddist
&& ptr
[0]); c
++)
241 if (sscanf(buf
, "%lf", &x1
) != 1)
253 atoms
->nres
= newres
+ 1;
256 fgets2(line
, STRLEN
, in
);
257 if (sscanf(line
, "%lf%lf%lf", &x1
, &y1
, &z1
) != 3)
259 gmx_warning("Bad box in file %s", infile
);
261 /* Generate a cubic box */
262 for (m
= 0; (m
< DIM
); m
++)
264 xmin
[m
] = xmax
[m
] = x
[0][m
];
266 for (i
= 1; (i
< atoms
->nr
); i
++)
268 for (m
= 0; (m
< DIM
); m
++)
270 xmin
[m
] = std::min(xmin
[m
], x
[i
][m
]);
271 xmax
[m
] = std::max(xmax
[m
], x
[i
][m
]);
274 for (i
= 0; i
< DIM
; i
++)
276 for (m
= 0; m
< DIM
; m
++)
281 for (m
= 0; (m
< DIM
); m
++)
283 box
[m
][m
] = (xmax
[m
]-xmin
[m
]);
285 fprintf(stderr
, "Generated a cubic box %8.3f x %8.3f x %8.3f\n",
286 box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
]);
290 /* We found the first three values, the diagonal elements */
294 if (sscanf (line
, "%*f%*f%*f%lf%lf%lf%lf%lf%lf",
295 &x1
, &y1
, &z1
, &x2
, &y2
, &z2
) != 6)
297 x1
= y1
= z1
= x2
= y2
= z2
= 0.0;
310 void gmx_gro_read_conf(const char *infile
,
311 t_symtab
*symtab
, char ***name
, t_atoms
*atoms
,
312 rvec x
[], rvec
*v
, matrix box
)
314 FILE *in
= gmx_fio_fopen(infile
, "r");
317 get_w_conf(in
, infile
, title
, symtab
, atoms
, &ndec
, x
, v
, box
);
318 *name
= put_symtab(symtab
, title
);
322 static gmx_bool
gmx_one_before_eof(FILE *fp
)
327 if ((beof
= fread(data
, 1, 1, fp
)) == 1)
329 gmx_fseek(fp
, -1, SEEK_CUR
);
334 gmx_bool
gro_next_x_or_v(FILE *status
, t_trxframe
*fr
)
338 char title
[STRLEN
], *p
;
342 if (gmx_one_before_eof(status
))
347 open_symtab(&symtab
);
348 atoms
.nr
= fr
->natoms
;
349 snew(atoms
.atom
, fr
->natoms
);
350 atoms
.nres
= fr
->natoms
;
351 snew(atoms
.resinfo
, fr
->natoms
);
352 snew(atoms
.atomname
, fr
->natoms
);
354 fr
->bV
= get_w_conf(status
, title
, title
, &symtab
, &atoms
, &ndec
, fr
->x
, fr
->v
, fr
->box
);
357 /* prec = 10^ndec: */
358 for (i
= 0; i
< ndec
; i
++)
368 sfree(atoms
.resinfo
);
369 sfree(atoms
.atomname
);
370 done_symtab(&symtab
);
372 if ((p
= strstr(title
, "t=")) != NULL
)
375 if (sscanf(p
, "%lf", &tt
) == 1)
387 if (atoms
.nr
!= fr
->natoms
)
389 gmx_fatal(FARGS
, "Number of atoms in gro frame (%d) doesn't match the number in the previous frame (%d)", atoms
.nr
, fr
->natoms
);
395 int gro_first_x_or_v(FILE *status
, t_trxframe
*fr
)
400 fprintf(stderr
, "Reading frames from gro file");
401 get_coordnum_fp(status
, title
, &fr
->natoms
);
403 fprintf(stderr
, " '%s', %d atoms.\n", title
, fr
->natoms
);
408 gmx_file("No coordinates in gro file");
411 snew(fr
->x
, fr
->natoms
);
412 snew(fr
->v
, fr
->natoms
);
413 gro_next_x_or_v(status
, fr
);
418 static const char *get_hconf_format(bool haveVelocities
)
422 return "%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f\n";
426 return "%8.3f%8.3f%8.3f\n";
431 static void write_hconf_box(FILE *out
, const matrix box
)
433 if (box
[XX
][YY
] || box
[XX
][ZZ
] || box
[YY
][XX
] || box
[YY
][ZZ
] ||
434 box
[ZZ
][XX
] || box
[ZZ
][YY
])
436 fprintf(out
, "%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n",
437 box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
],
438 box
[XX
][YY
], box
[XX
][ZZ
], box
[YY
][XX
],
439 box
[YY
][ZZ
], box
[ZZ
][XX
], box
[ZZ
][YY
]);
443 fprintf(out
, "%10.5f%10.5f%10.5f\n",
444 box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
]);
448 void write_hconf_indexed_p(FILE *out
, const char *title
, const t_atoms
*atoms
,
449 int nx
, const int index
[],
450 const rvec
*x
, const rvec
*v
, const matrix box
)
452 char resnm
[6], nm
[6];
453 int ai
, i
, resind
, resnr
;
455 fprintf(out
, "%s\n", (title
&& title
[0]) ? title
: gmx::bromacs().c_str());
456 fprintf(out
, "%5d\n", nx
);
458 const char *format
= get_hconf_format(v
!= NULL
);
460 for (i
= 0; (i
< nx
); i
++)
464 resind
= atoms
->atom
[ai
].resind
;
465 std::strncpy(resnm
, " ??? ", sizeof(resnm
)-1);
466 if (resind
< atoms
->nres
)
468 std::strncpy(resnm
, *atoms
->resinfo
[resind
].name
, sizeof(resnm
)-1);
469 resnr
= atoms
->resinfo
[resind
].nr
;
473 std::strncpy(resnm
, " ??? ", sizeof(resnm
)-1);
479 std::strncpy(nm
, *atoms
->atomname
[ai
], sizeof(nm
)-1);
483 std::strncpy(nm
, " ??? ", sizeof(nm
)-1);
486 fprintf(out
, "%5d%-5.5s%5.5s%5d", resnr
%100000, resnm
, nm
, (ai
+1)%100000);
487 /* next fprintf uses built format string */
491 x
[ai
][XX
], x
[ai
][YY
], x
[ai
][ZZ
], v
[ai
][XX
], v
[ai
][YY
], v
[ai
][ZZ
]);
496 x
[ai
][XX
], x
[ai
][YY
], x
[ai
][ZZ
]);
500 write_hconf_box(out
, box
);
505 void write_hconf_mtop(FILE *out
, const char *title
, gmx_mtop_t
*mtop
,
506 const rvec
*x
, const rvec
*v
, const matrix box
)
509 gmx_mtop_atomloop_all_t aloop
;
511 char *atomname
, *resname
;
513 fprintf(out
, "%s\n", (title
&& title
[0]) ? title
: gmx::bromacs().c_str());
514 fprintf(out
, "%5d\n", mtop
->natoms
);
516 const char *format
= get_hconf_format(v
!= NULL
);
518 aloop
= gmx_mtop_atomloop_all_init(mtop
);
519 while (gmx_mtop_atomloop_all_next(aloop
, &i
, &atom
))
521 gmx_mtop_atomloop_all_names(aloop
, &atomname
, &resnr
, &resname
);
523 fprintf(out
, "%5d%-5.5s%5.5s%5d",
524 resnr
%100000, resname
, atomname
, (i
+1)%100000);
525 /* next fprintf uses built format string */
529 x
[i
][XX
], x
[i
][YY
], x
[i
][ZZ
], v
[i
][XX
], v
[i
][YY
], v
[i
][ZZ
]);
534 x
[i
][XX
], x
[i
][YY
], x
[i
][ZZ
]);
538 write_hconf_box(out
, box
);
543 void write_hconf_p(FILE *out
, const char *title
, const t_atoms
*atoms
,
544 const rvec
*x
, const rvec
*v
, const matrix box
)
550 for (i
= 0; (i
< atoms
->nr
); i
++)
554 write_hconf_indexed_p(out
, title
, atoms
, atoms
->nr
, aa
, x
, v
, box
);
558 void write_conf_p(const char *outfile
, const char *title
,
559 const t_atoms
*atoms
,
560 const rvec
*x
, const rvec
*v
, const matrix box
)
564 out
= gmx_fio_fopen(outfile
, "w");
565 write_hconf_p(out
, title
, atoms
, x
, v
, box
);