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.
44 #include "gromacs/fileio/trxio.h"
45 #include "gromacs/math/vec.h"
46 #include "gromacs/topology/atoms.h"
47 #include "gromacs/topology/symtab.h"
48 #include "gromacs/trajectory/trajectoryframe.h"
49 #include "gromacs/utility/cstringutil.h"
50 #include "gromacs/utility/fatalerror.h"
51 #include "gromacs/utility/gmxassert.h"
52 #include "gromacs/utility/smalloc.h"
56 static int read_g96_pos(char line
[], t_symtab
*symtab
,
57 FILE *fp
, const char *infile
,
62 int nwanted
, natoms
, atnr
, resnr
= 0, oldres
, newres
, shift
;
63 char anm
[STRLEN
], resnm
[STRLEN
];
69 if (fr
->atoms
!= nullptr)
71 GMX_RELEASE_ASSERT(symtab
!= nullptr, "Reading a conformation from a g96 format with atom data requires a valid symbol table");
76 atoms
->haveMass
= FALSE
;
77 atoms
->haveCharge
= FALSE
;
78 atoms
->haveType
= FALSE
;
79 atoms
->haveBState
= FALSE
;
80 atoms
->havePdbInfo
= FALSE
;
96 oldres
= -666; /* Unlikely number for the first residue! */
98 while (!bEnd
&& fgets2(line
, STRLEN
, fp
))
100 bEnd
= (std::strncmp(line
, "END", 3) == 0);
101 if (!bEnd
&& (line
[0] != '#'))
103 if (sscanf(line
+shift
, "%15lf%15lf%15lf", &db1
, &db2
, &db3
) != 3)
105 gmx_fatal(FARGS
, "Did not find 3 coordinates for atom %d in %s\n",
108 if ((nwanted
!= -1) && (natoms
>= nwanted
))
111 "Found more coordinates (%d) in %s than expected %d\n",
112 natoms
, infile
, nwanted
);
117 (sscanf(line
, "%5d%c%5s%c%5s%7d", &resnr
, &c1
, resnm
, &c2
, anm
, &atnr
)
127 strncpy(resnm
, "???", sizeof(resnm
)-1);
129 strncpy(anm
, "???", sizeof(anm
)-1);
131 atoms
->atomname
[natoms
] = put_symtab(symtab
, anm
);
136 if (newres
>= atoms
->nr
)
138 gmx_fatal(FARGS
, "More residues than atoms in %s (natoms = %d)",
141 atoms
->atom
[natoms
].resind
= newres
;
142 if (newres
+1 > atoms
->nres
)
144 atoms
->nres
= newres
+1;
146 t_atoms_set_resinfo(atoms
, natoms
, symtab
, resnm
, resnr
, ' ', 0, ' ');
150 atoms
->atom
[natoms
].resind
= newres
;
155 fr
->x
[natoms
][0] = db1
;
156 fr
->x
[natoms
][1] = db2
;
157 fr
->x
[natoms
][2] = db3
;
162 if ((nwanted
!= -1) && natoms
!= nwanted
)
165 "Warning: found less coordinates (%d) in %s than expected %d\n",
166 natoms
, infile
, nwanted
);
175 static int read_g96_vel(char line
[], FILE *fp
, const char *infile
,
179 int nwanted
, natoms
= -1, shift
;
180 double db1
, db2
, db3
;
182 nwanted
= fr
->natoms
;
186 if (strcmp(line
, "VELOCITYRED") == 0)
196 while (!bEnd
&& fgets2(line
, STRLEN
, fp
))
198 bEnd
= (strncmp(line
, "END", 3) == 0);
199 if (!bEnd
&& (line
[0] != '#'))
201 if (sscanf(line
+shift
, "%15lf%15lf%15lf", &db1
, &db2
, &db3
) != 3)
203 gmx_fatal(FARGS
, "Did not find 3 velocities for atom %d in %s",
206 if ((nwanted
!= -1) && (natoms
>= nwanted
))
208 gmx_fatal(FARGS
, "Found more velocities (%d) in %s than expected %d\n",
209 natoms
, infile
, nwanted
);
213 fr
->v
[natoms
][0] = db1
;
214 fr
->v
[natoms
][1] = db2
;
215 fr
->v
[natoms
][2] = db3
;
220 if ((nwanted
!= -1) && (natoms
!= nwanted
))
223 "Warning: found less velocities (%d) in %s than expected %d\n",
224 natoms
, infile
, nwanted
);
231 int read_g96_conf(FILE *fp
, const char *infile
, char **name
, t_trxframe
*fr
,
232 t_symtab
*symtab
, char *line
)
234 gmx_bool bAtStart
, bTime
, bAtoms
, bPos
, bVel
, bBox
, bEnd
, bFinished
;
236 double db1
, db2
, db3
, db4
, db5
, db6
, db7
, db8
, db9
;
238 bAtStart
= (ftell(fp
) == 0);
240 clear_trxframe(fr
, FALSE
);
246 bool foundTitle
= false;
247 while (!foundTitle
&& fgets2(line
, STRLEN
, fp
))
249 foundTitle
= (std::strcmp(line
, "TITLE") == 0);
251 fgets2(line
, STRLEN
, fp
);
254 *name
= gmx_strdup(line
);
257 while (!bEnd
&& fgets2(line
, STRLEN
, fp
))
259 bEnd
= (std::strcmp(line
, "END") == 0);
261 fgets2(line
, STRLEN
, fp
);
264 /* Do not get a line if we are not at the start of the file, *
265 * because without a parameter file we don't know what is in *
266 * the trajectory and we have already read the line in the *
267 * previous call (VERY DIRTY). */
271 bTime
= (std::strcmp(line
, "TIMESTEP") == 0);
272 bAtoms
= (std::strcmp(line
, "POSITION") == 0);
273 bPos
= (bAtoms
|| (strcmp(line
, "POSITIONRED") == 0));
274 bVel
= (std::strncmp(line
, "VELOCITY", 8) == 0);
275 bBox
= (std::strcmp(line
, "BOX") == 0);
278 if (!fr
->bTime
&& !fr
->bX
)
284 bFinished
= (fgets2(line
, STRLEN
, fp
) == nullptr);
286 while (!bFinished
&& (line
[0] == '#'));
287 sscanf(line
, "%15" SCNd64
"%15lf", &(fr
->step
), &db1
);
301 natoms
= read_g96_pos(line
, symtab
, fp
, infile
, fr
);
311 natoms
= read_g96_vel(line
, fp
, infile
, fr
);
318 while (!bEnd
&& fgets2(line
, STRLEN
, fp
))
320 bEnd
= (strncmp(line
, "END", 3) == 0);
321 if (!bEnd
&& (line
[0] != '#'))
323 nbp
= sscanf(line
, "%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf",
324 &db1
, &db2
, &db3
, &db4
, &db5
, &db6
, &db7
, &db8
, &db9
);
327 gmx_fatal(FARGS
, "Found a BOX line, but no box in %s", infile
);
329 fr
->box
[XX
][XX
] = db1
;
330 fr
->box
[YY
][YY
] = db2
;
331 fr
->box
[ZZ
][ZZ
] = db3
;
334 fr
->box
[XX
][YY
] = db4
;
335 fr
->box
[XX
][ZZ
] = db5
;
336 fr
->box
[YY
][XX
] = db6
;
337 fr
->box
[YY
][ZZ
] = db7
;
338 fr
->box
[ZZ
][XX
] = db8
;
339 fr
->box
[ZZ
][YY
] = db9
;
346 while (!bFinished
&& (fgets2(line
, STRLEN
, fp
) != nullptr));
353 void write_g96_conf(FILE *out
, const char *title
, const t_trxframe
*fr
,
354 int nindex
, const int *index
)
370 fprintf(out
, "TITLE\n%s\nEND\n", title
);
371 if (fr
->bStep
|| fr
->bTime
)
373 /* Officially the time format is %15.9, which is not enough for 10 ns */
374 fprintf(out
, "TIMESTEP\n%15" PRId64
"%15.6f\nEND\n", fr
->step
, fr
->time
);
380 fprintf(out
, "POSITION\n");
381 for (i
= 0; i
< nout
; i
++)
391 fprintf(out
, "%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
392 (atoms
->resinfo
[atoms
->atom
[a
].resind
].nr
) % 100000,
393 *atoms
->resinfo
[atoms
->atom
[a
].resind
].name
,
394 *atoms
->atomname
[a
], (i
+1) % 10000000,
395 fr
->x
[a
][XX
], fr
->x
[a
][YY
], fr
->x
[a
][ZZ
]);
400 fprintf(out
, "POSITIONRED\n");
401 for (i
= 0; i
< nout
; i
++)
411 fprintf(out
, "%15.9f%15.9f%15.9f\n",
412 fr
->x
[a
][XX
], fr
->x
[a
][YY
], fr
->x
[a
][ZZ
]);
415 fprintf(out
, "END\n");
421 fprintf(out
, "VELOCITY\n");
422 for (i
= 0; i
< nout
; i
++)
432 fprintf(out
, "%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
433 (atoms
->resinfo
[atoms
->atom
[a
].resind
].nr
) % 100000,
434 *atoms
->resinfo
[atoms
->atom
[a
].resind
].name
,
435 *atoms
->atomname
[a
], (i
+1) % 10000000,
436 fr
->v
[a
][XX
], fr
->v
[a
][YY
], fr
->v
[a
][ZZ
]);
441 fprintf(out
, "VELOCITYRED\n");
442 for (i
= 0; i
< nout
; i
++)
452 fprintf(out
, "%15.9f%15.9f%15.9f\n",
453 fr
->v
[a
][XX
], fr
->v
[a
][YY
], fr
->v
[a
][ZZ
]);
456 fprintf(out
, "END\n");
460 fprintf(out
, "BOX\n");
461 fprintf(out
, "%15.9f%15.9f%15.9f",
462 fr
->box
[XX
][XX
], fr
->box
[YY
][YY
], fr
->box
[ZZ
][ZZ
]);
463 if ((fr
->box
[XX
][YY
] != 0.0F
) || (fr
->box
[XX
][ZZ
] != 0.0F
) || (fr
->box
[YY
][XX
] != 0.0F
) ||
464 (fr
->box
[YY
][ZZ
] != 0.0F
) || (fr
->box
[ZZ
][XX
] != 0.0F
) || (fr
->box
[ZZ
][YY
] != 0.0F
))
466 fprintf(out
, "%15.9f%15.9f%15.9f%15.9f%15.9f%15.9f",
467 fr
->box
[XX
][YY
], fr
->box
[XX
][ZZ
], fr
->box
[YY
][XX
],
468 fr
->box
[YY
][ZZ
], fr
->box
[ZZ
][XX
], fr
->box
[ZZ
][YY
]);
471 fprintf(out
, "END\n");