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 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.
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.
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.
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.
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.
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.
48 #include "gromacs/fileio/gmxfio.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/trajectory/trajectoryframe.h"
54 #include "gromacs/utility/coolstuff.h"
55 #include "gromacs/utility/cstringutil.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/futil.h"
58 #include "gromacs/utility/smalloc.h"
60 static void get_coordnum_fp(FILE* in
, char* title
, int* natoms
)
62 char line
[STRLEN
+ 1];
64 fgets2(title
, STRLEN
, in
);
65 fgets2(line
, STRLEN
, in
);
66 if (sscanf(line
, "%d", natoms
) != 1)
68 gmx_fatal(FARGS
, "gro file does not have the number of atoms on the second line");
72 void get_coordnum(const char* infile
, int* natoms
)
77 in
= gmx_fio_fopen(infile
, "r");
78 get_coordnum_fp(in
, title
, natoms
);
82 /* Note that the .gro reading routine still support variable precision
83 * for backward compatibility with old .gro files.
84 * We have removed writing of variable precision to avoid compatibility
85 * issues with other software packages.
87 static gmx_bool
get_w_conf(FILE* in
,
98 char resname
[6], oldresname
[6];
99 char line
[STRLEN
+ 1], *ptr
;
101 double x1
, y1
, z1
, x2
, y2
, z2
;
103 int natoms
, i
, m
, resnr
, newres
, oldres
, ddist
, c
;
104 gmx_bool bFirst
, bVel
, oldResFirst
;
112 /* Read the title and number of atoms */
113 get_coordnum_fp(in
, title
, &natoms
);
115 if (natoms
> atoms
->nr
)
117 gmx_fatal(FARGS
, "gro file contains more atoms (%d) than expected (%d)", natoms
, atoms
->nr
);
119 else if (natoms
< atoms
->nr
)
122 "Warning: gro file contains less atoms (%d) than expected"
127 atoms
->haveMass
= FALSE
;
128 atoms
->haveCharge
= FALSE
;
129 atoms
->haveType
= FALSE
;
130 atoms
->haveBState
= FALSE
;
131 atoms
->havePdbInfo
= FALSE
;
138 oldresname
[0] = '\0';
140 /* just pray the arrays are big enough */
141 for (i
= 0; (i
< natoms
); i
++)
143 if ((fgets2(line
, STRLEN
, in
)) == nullptr)
145 gmx_fatal(FARGS
, "Unexpected end of file in file %s at line %d", infile
, i
+ 2);
147 if (strlen(line
) < 39)
149 gmx_fatal(FARGS
, "Invalid line in %s for atom %d:\n%s", infile
, i
+ 1, line
);
152 /* determine read precision from distance between periods
157 p1
= strchr(line
, '.');
160 gmx_fatal(FARGS
, "A coordinate in file %s does not contain a '.'", infile
);
162 p2
= strchr(&p1
[1], '.');
165 gmx_fatal(FARGS
, "A coordinate in file %s does not contain a '.'", infile
);
170 p3
= strchr(&p2
[1], '.');
173 gmx_fatal(FARGS
, "A coordinate in file %s does not contain a '.'", infile
);
176 if (p3
- p2
!= ddist
)
179 "The spacing of the decimal points in file %s is not consistent for x, y "
186 memcpy(name
, line
, 5);
188 sscanf(name
, "%d", &resnr
);
189 sscanf(line
+ 5, "%5s", resname
);
191 if (!oldResFirst
|| oldres
!= resnr
|| strncmp(resname
, oldresname
, sizeof(resname
)) != 0)
196 if (newres
>= natoms
)
198 gmx_fatal(FARGS
, "More residues than atoms in %s (natoms = %d)", infile
, natoms
);
200 atoms
->atom
[i
].resind
= newres
;
201 t_atoms_set_resinfo(atoms
, i
, symtab
, resname
, resnr
, ' ', 0, ' ');
205 atoms
->atom
[i
].resind
= newres
;
209 std::memcpy(name
, line
+ 10, 5);
210 atoms
->atomname
[i
] = put_symtab(symtab
, name
);
212 /* Copy resname to oldresname after we are done with the sanity check above */
213 std::strncpy(oldresname
, resname
, sizeof(oldresname
));
215 /* eventueel controle atomnumber met i+1 */
217 /* coordinates (start after residue data) */
219 /* Read fixed format */
220 for (m
= 0; m
< DIM
; m
++)
222 for (c
= 0; (c
< ddist
&& ptr
[0]); c
++)
228 if (sscanf(buf
, "%lf %lf", &x1
, &x2
) != 1)
231 "Something is wrong in the coordinate formatting of file %s. Note that "
232 "gro is fixed format (see the manual)",
241 /* velocities (start after residues and coordinates) */
244 /* Read fixed format */
245 for (m
= 0; m
< DIM
; m
++)
247 for (c
= 0; (c
< ddist
&& ptr
[0]); c
++)
253 if (sscanf(buf
, "%lf", &x1
) != 1)
265 atoms
->nres
= newres
+ 1;
268 fgets2(line
, STRLEN
, in
);
269 if (sscanf(line
, "%lf%lf%lf", &x1
, &y1
, &z1
) != 3)
271 gmx_warning("Bad box in file %s", infile
);
273 /* Generate a cubic box */
274 for (m
= 0; (m
< DIM
); m
++)
276 xmin
[m
] = xmax
[m
] = x
[0][m
];
278 for (i
= 1; (i
< atoms
->nr
); i
++)
280 for (m
= 0; (m
< DIM
); m
++)
282 xmin
[m
] = std::min(xmin
[m
], x
[i
][m
]);
283 xmax
[m
] = std::max(xmax
[m
], x
[i
][m
]);
286 for (i
= 0; i
< DIM
; i
++)
288 for (m
= 0; m
< DIM
; m
++)
293 for (m
= 0; (m
< DIM
); m
++)
295 box
[m
][m
] = (xmax
[m
] - xmin
[m
]);
297 fprintf(stderr
, "Generated a cubic box %8.3f x %8.3f x %8.3f\n", box
[XX
][XX
], box
[YY
][YY
],
302 /* We found the first three values, the diagonal elements */
306 if (sscanf(line
, "%*f%*f%*f%lf%lf%lf%lf%lf%lf", &x1
, &y1
, &z1
, &x2
, &y2
, &z2
) != 6)
308 x1
= y1
= z1
= x2
= y2
= z2
= 0.0;
321 void gmx_gro_read_conf(const char* infile
, t_symtab
* symtab
, char** name
, t_atoms
* atoms
, rvec x
[], rvec
* v
, matrix box
)
323 FILE* in
= gmx_fio_fopen(infile
, "r");
326 get_w_conf(in
, infile
, title
, symtab
, atoms
, &ndec
, x
, v
, box
);
329 *name
= gmx_strdup(title
);
334 static gmx_bool
gmx_one_before_eof(FILE* fp
)
337 gmx_bool beof
= fread(data
, 1, 1, fp
) != 1;
341 gmx_fseek(fp
, -1, SEEK_CUR
);
346 gmx_bool
gro_next_x_or_v(FILE* status
, t_trxframe
* fr
)
350 char title
[STRLEN
], *p
;
354 if (gmx_one_before_eof(status
))
359 open_symtab(&symtab
);
360 atoms
.nr
= fr
->natoms
;
361 snew(atoms
.atom
, fr
->natoms
);
362 atoms
.nres
= fr
->natoms
;
363 snew(atoms
.resinfo
, fr
->natoms
);
364 snew(atoms
.atomname
, fr
->natoms
);
366 fr
->bV
= get_w_conf(status
, title
, title
, &symtab
, &atoms
, &ndec
, fr
->x
, fr
->v
, fr
->box
);
369 /* prec = 10^ndec: */
370 for (i
= 0; i
< ndec
; i
++)
378 sfree(atoms
.resinfo
);
379 sfree(atoms
.atomname
);
380 done_symtab(&symtab
);
382 if ((p
= strstr(title
, "t=")) != nullptr)
385 if (sscanf(p
, "%lf", &tt
) == 1)
397 if ((p
= std::strstr(title
, "step=")) != nullptr)
400 fr
->step
= 0; // Default value if fr-bStep is false
401 fr
->bStep
= (sscanf(p
, "%" SCNd64
, &fr
->step
) == 1);
404 if (atoms
.nr
!= fr
->natoms
)
407 "Number of atoms in gro frame (%d) doesn't match the number in the previous "
409 atoms
.nr
, fr
->natoms
);
415 int gro_first_x_or_v(FILE* status
, t_trxframe
* fr
)
420 fprintf(stderr
, "Reading frames from gro file");
421 get_coordnum_fp(status
, title
, &fr
->natoms
);
423 fprintf(stderr
, " '%s', %d atoms.\n", title
, fr
->natoms
);
426 gmx_file("No coordinates in gro file");
429 snew(fr
->x
, fr
->natoms
);
430 snew(fr
->v
, fr
->natoms
);
431 gro_next_x_or_v(status
, fr
);
436 static const char* get_hconf_format(bool haveVelocities
)
440 return "%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f\n";
444 return "%8.3f%8.3f%8.3f\n";
448 static void write_hconf_box(FILE* out
, const matrix box
)
450 if ((box
[XX
][YY
] != 0.0F
) || (box
[XX
][ZZ
] != 0.0F
) || (box
[YY
][XX
] != 0.0F
)
451 || (box
[YY
][ZZ
] != 0.0F
) || (box
[ZZ
][XX
] != 0.0F
) || (box
[ZZ
][YY
] != 0.0F
))
453 fprintf(out
, "%10.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n", box
[XX
][XX
],
454 box
[YY
][YY
], box
[ZZ
][ZZ
], box
[XX
][YY
], box
[XX
][ZZ
], box
[YY
][XX
], box
[YY
][ZZ
],
455 box
[ZZ
][XX
], box
[ZZ
][YY
]);
459 fprintf(out
, "%10.5f %9.5f %9.5f\n", box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
]);
463 void write_hconf_indexed_p(FILE* out
,
465 const t_atoms
* atoms
,
472 int ai
, i
, resind
, resnr
;
474 fprintf(out
, "%s\n", (title
&& title
[0]) ? title
: gmx::bromacs().c_str());
475 fprintf(out
, "%5d\n", nx
);
477 const char* format
= get_hconf_format(v
!= nullptr);
479 for (i
= 0; (i
< nx
); i
++)
483 resind
= atoms
->atom
[ai
].resind
;
485 if (resind
< atoms
->nres
)
487 resnm
= *atoms
->resinfo
[resind
].name
;
488 resnr
= atoms
->resinfo
[resind
].nr
;
499 nm
= *atoms
->atomname
[ai
];
506 fprintf(out
, "%5d%-5.5s%5.5s%5d", resnr
% 100000, resnm
.c_str(), nm
.c_str(), (ai
+ 1) % 100000);
507 /* next fprintf uses built format string */
510 fprintf(out
, format
, x
[ai
][XX
], x
[ai
][YY
], x
[ai
][ZZ
], v
[ai
][XX
], v
[ai
][YY
], v
[ai
][ZZ
]);
514 fprintf(out
, format
, x
[ai
][XX
], x
[ai
][YY
], x
[ai
][ZZ
]);
518 write_hconf_box(out
, box
);
523 void write_hconf_mtop(FILE* out
, const char* title
, const gmx_mtop_t
* mtop
, const rvec
* x
, const rvec
* v
, const matrix box
)
525 fprintf(out
, "%s\n", (title
&& title
[0]) ? title
: gmx::bromacs().c_str());
526 fprintf(out
, "%5d\n", mtop
->natoms
);
528 const char* format
= get_hconf_format(v
!= nullptr);
530 for (const AtomProxy atomP
: AtomRange(*mtop
))
532 int i
= atomP
.globalAtomNumber();
533 int residueNumber
= atomP
.residueNumber();
534 const char* atomName
= atomP
.atomName();
535 const char* residueName
= atomP
.residueName();
537 fprintf(out
, "%5d%-5.5s%5.5s%5d", residueNumber
% 100000, residueName
, atomName
, (i
+ 1) % 100000);
538 /* next fprintf uses built format string */
541 fprintf(out
, format
, x
[i
][XX
], x
[i
][YY
], x
[i
][ZZ
], v
[i
][XX
], v
[i
][YY
], v
[i
][ZZ
]);
545 fprintf(out
, format
, x
[i
][XX
], x
[i
][YY
], x
[i
][ZZ
]);
549 write_hconf_box(out
, box
);
554 void write_hconf_p(FILE* out
, const char* title
, const t_atoms
* atoms
, const rvec
* x
, const rvec
* v
, const matrix box
)
560 for (i
= 0; (i
< atoms
->nr
); i
++)
564 write_hconf_indexed_p(out
, title
, atoms
, atoms
->nr
, aa
, x
, v
, box
);
568 void write_conf_p(const char* outfile
,
570 const t_atoms
* atoms
,
577 out
= gmx_fio_fopen(outfile
, "w");
578 write_hconf_p(out
, title
, atoms
, x
, v
, box
);