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.
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
)) == nullptr)
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
)) != 0)
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
);
320 *name
= gmx_strdup(title
);
325 static gmx_bool
gmx_one_before_eof(FILE *fp
)
328 gmx_bool beof
= fread(data
, 1, 1, fp
) != 1;
332 gmx_fseek(fp
, -1, SEEK_CUR
);
337 gmx_bool
gro_next_x_or_v(FILE *status
, t_trxframe
*fr
)
341 char title
[STRLEN
], *p
;
345 if (gmx_one_before_eof(status
))
350 open_symtab(&symtab
);
351 atoms
.nr
= fr
->natoms
;
352 snew(atoms
.atom
, fr
->natoms
);
353 atoms
.nres
= fr
->natoms
;
354 snew(atoms
.resinfo
, fr
->natoms
);
355 snew(atoms
.atomname
, fr
->natoms
);
357 fr
->bV
= get_w_conf(status
, title
, title
, &symtab
, &atoms
, &ndec
, fr
->x
, fr
->v
, fr
->box
);
360 /* prec = 10^ndec: */
361 for (i
= 0; i
< ndec
; i
++)
369 sfree(atoms
.resinfo
);
370 sfree(atoms
.atomname
);
371 done_symtab(&symtab
);
373 if ((p
= strstr(title
, "t=")) != nullptr)
376 if (sscanf(p
, "%lf", &tt
) == 1)
388 if ((p
= std::strstr(title
, "step=")) != nullptr)
391 fr
->step
= 0; // Default value if fr-bStep is false
392 fr
->bStep
= (sscanf(p
, "%" SCNd64
, &fr
->step
) == 1);
395 if (atoms
.nr
!= fr
->natoms
)
397 gmx_fatal(FARGS
, "Number of atoms in gro frame (%d) doesn't match the number in the previous frame (%d)", atoms
.nr
, fr
->natoms
);
403 int gro_first_x_or_v(FILE *status
, t_trxframe
*fr
)
408 fprintf(stderr
, "Reading frames from gro file");
409 get_coordnum_fp(status
, title
, &fr
->natoms
);
411 fprintf(stderr
, " '%s', %d atoms.\n", title
, fr
->natoms
);
414 gmx_file("No coordinates in gro file");
417 snew(fr
->x
, fr
->natoms
);
418 snew(fr
->v
, fr
->natoms
);
419 gro_next_x_or_v(status
, fr
);
424 static const char *get_hconf_format(bool haveVelocities
)
428 return "%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f\n";
432 return "%8.3f%8.3f%8.3f\n";
437 static void write_hconf_box(FILE *out
, const matrix box
)
439 if ((box
[XX
][YY
] != 0.0F
) || (box
[XX
][ZZ
] != 0.0F
) || (box
[YY
][XX
] != 0.0F
) || (box
[YY
][ZZ
] != 0.0F
) ||
440 (box
[ZZ
][XX
] != 0.0F
) || (box
[ZZ
][YY
] != 0.0F
))
442 fprintf(out
, "%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n",
443 box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
],
444 box
[XX
][YY
], box
[XX
][ZZ
], box
[YY
][XX
],
445 box
[YY
][ZZ
], box
[ZZ
][XX
], box
[ZZ
][YY
]);
449 fprintf(out
, "%10.5f%10.5f%10.5f\n",
450 box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
]);
454 void write_hconf_indexed_p(FILE *out
, const char *title
, const t_atoms
*atoms
,
455 int nx
, const int index
[],
456 const rvec
*x
, const rvec
*v
, const matrix box
)
458 int ai
, i
, resind
, resnr
;
460 fprintf(out
, "%s\n", (title
&& title
[0]) ? title
: gmx::bromacs().c_str());
461 fprintf(out
, "%5d\n", nx
);
463 const char *format
= get_hconf_format(v
!= nullptr);
465 for (i
= 0; (i
< nx
); i
++)
469 resind
= atoms
->atom
[ai
].resind
;
471 if (resind
< atoms
->nres
)
473 resnm
= *atoms
->resinfo
[resind
].name
;
474 resnr
= atoms
->resinfo
[resind
].nr
;
485 nm
= *atoms
->atomname
[i
];
492 fprintf(out
, "%5d%-5.5s%5.5s%5d", resnr
%100000, resnm
.c_str(), nm
.c_str(), (ai
+1)%100000);
493 /* next fprintf uses built format string */
497 x
[ai
][XX
], x
[ai
][YY
], x
[ai
][ZZ
], v
[ai
][XX
], v
[ai
][YY
], v
[ai
][ZZ
]);
502 x
[ai
][XX
], x
[ai
][YY
], x
[ai
][ZZ
]);
506 write_hconf_box(out
, box
);
511 void write_hconf_mtop(FILE *out
, const char *title
, const gmx_mtop_t
*mtop
,
512 const rvec
*x
, const rvec
*v
, const matrix box
)
514 fprintf(out
, "%s\n", (title
&& title
[0]) ? title
: gmx::bromacs().c_str());
515 fprintf(out
, "%5d\n", mtop
->natoms
);
517 const char *format
= get_hconf_format(v
!= nullptr);
519 for (const AtomProxy atomP
: AtomRange(*mtop
))
521 int i
= atomP
.globalAtomNumber();
522 int residueNumber
= atomP
.residueNumber();
523 const char *atomName
= atomP
.atomName();
524 const char *residueName
= atomP
.residueName();
526 fprintf(out
, "%5d%-5.5s%5.5s%5d",
527 residueNumber
%100000, residueName
, atomName
, (i
+1)%100000);
528 /* next fprintf uses built format string */
532 x
[i
][XX
], x
[i
][YY
], x
[i
][ZZ
], v
[i
][XX
], v
[i
][YY
], v
[i
][ZZ
]);
537 x
[i
][XX
], x
[i
][YY
], x
[i
][ZZ
]);
541 write_hconf_box(out
, box
);
546 void write_hconf_p(FILE *out
, const char *title
, const t_atoms
*atoms
,
547 const rvec
*x
, const rvec
*v
, const matrix box
)
553 for (i
= 0; (i
< atoms
->nr
); i
++)
557 write_hconf_indexed_p(out
, title
, atoms
, atoms
->nr
, aa
, x
, v
, box
);
561 void write_conf_p(const char *outfile
, const char *title
,
562 const t_atoms
*atoms
,
563 const rvec
*x
, const rvec
*v
, const matrix box
)
567 out
= gmx_fio_fopen(outfile
, "w");
568 write_hconf_p(out
, title
, atoms
, x
, v
, box
);