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.
44 #include "gromacs/fileio/espio.h"
45 #include "gromacs/fileio/filetypes.h"
46 #include "gromacs/fileio/g96io.h"
47 #include "gromacs/fileio/gmxfio.h"
48 #include "gromacs/fileio/groio.h"
49 #include "gromacs/fileio/pdbio.h"
50 #include "gromacs/fileio/tpxio.h"
51 #include "gromacs/fileio/trxio.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/topology/atoms.h"
54 #include "gromacs/topology/block.h"
55 #include "gromacs/topology/mtop_util.h"
56 #include "gromacs/topology/symtab.h"
57 #include "gromacs/topology/topology.h"
58 #include "gromacs/trajectory/trajectoryframe.h"
59 #include "gromacs/utility/cstringutil.h"
60 #include "gromacs/utility/fatalerror.h"
61 #include "gromacs/utility/gmxassert.h"
62 #include "gromacs/utility/smalloc.h"
64 void write_sto_conf_indexed(const char *outfile
, const char *title
,
66 const rvec x
[], const rvec
*v
, int ePBC
, const matrix box
,
67 int nindex
, int index
[])
73 ftp
= fn2ftp(outfile
);
77 out
= gmx_fio_fopen(outfile
, "w");
78 write_hconf_indexed_p(out
, title
, atoms
, nindex
, index
, x
, v
, box
);
82 clear_trxframe(&fr
, TRUE
);
85 fr
.natoms
= atoms
->nr
;
87 fr
.atoms
= const_cast<t_atoms
*>(atoms
);
89 fr
.x
= const_cast<rvec
*>(x
);
93 fr
.v
= const_cast<rvec
*>(v
);
96 copy_mat(box
, fr
.box
);
97 out
= gmx_fio_fopen(outfile
, "w");
98 write_g96_conf(out
, &fr
, nindex
, index
);
105 out
= gmx_fio_fopen(outfile
, "w");
106 write_pdbfile_indexed(out
, title
, atoms
, x
, ePBC
, box
, ' ', -1, nindex
, index
, NULL
, TRUE
);
110 out
= gmx_fio_fopen(outfile
, "w");
111 write_espresso_conf_indexed(out
, title
, atoms
, nindex
, index
, x
, v
, box
);
115 gmx_fatal(FARGS
, "Sorry, can not write a topology to %s", outfile
);
118 gmx_incons("Not supported in write_sto_conf_indexed");
122 void write_sto_conf(const char *outfile
, const char *title
, const t_atoms
*atoms
,
123 const rvec x
[], const rvec
*v
, int ePBC
, const matrix box
)
129 ftp
= fn2ftp(outfile
);
133 write_conf_p(outfile
, title
, atoms
, x
, v
, box
);
136 clear_trxframe(&fr
, TRUE
);
139 fr
.natoms
= atoms
->nr
;
141 fr
.atoms
= const_cast<t_atoms
*>(atoms
); // TODO check
143 fr
.x
= const_cast<rvec
*>(x
);
147 fr
.v
= const_cast<rvec
*>(v
);
150 copy_mat(box
, fr
.box
);
151 out
= gmx_fio_fopen(outfile
, "w");
152 write_g96_conf(out
, &fr
, -1, NULL
);
158 out
= gmx_fio_fopen(outfile
, "w");
159 write_pdbfile(out
, title
, atoms
, x
, ePBC
, box
, ' ', -1, NULL
, TRUE
);
163 out
= gmx_fio_fopen(outfile
, "w");
164 write_espresso_conf_indexed(out
, title
, atoms
, atoms
->nr
, NULL
, x
, v
, box
);
168 gmx_fatal(FARGS
, "Sorry, can not write a topology to %s", outfile
);
171 gmx_incons("Not supported in write_sto_conf");
175 void write_sto_conf_mtop(const char *outfile
, const char *title
,
177 const rvec x
[], const rvec
*v
, int ePBC
, const matrix box
)
183 ftp
= fn2ftp(outfile
);
187 out
= gmx_fio_fopen(outfile
, "w");
188 write_hconf_mtop(out
, title
, mtop
, x
, v
, box
);
192 /* This is a brute force approach which requires a lot of memory.
193 * We should implement mtop versions of all writing routines.
195 atoms
= gmx_mtop_global_atoms(mtop
);
197 write_sto_conf(outfile
, title
, &atoms
, x
, v
, ePBC
, box
);
204 static void get_stx_coordnum(const char *infile
, int *natoms
)
209 char g96_line
[STRLEN
+1];
211 ftp
= fn2ftp(infile
);
212 range_check(ftp
, 0, efNR
);
216 get_coordnum(infile
, natoms
);
220 in
= gmx_fio_fopen(infile
, "r");
227 *natoms
= read_g96_conf(in
, infile
, &fr
, NULL
, g96_line
);
228 sfree(const_cast<char *>(fr
.title
));
235 in
= gmx_fio_fopen(infile
, "r");
236 get_pdb_coordnum(in
, natoms
);
240 *natoms
= get_espresso_coordnum(infile
);
243 gmx_fatal(FARGS
, "File type %s not supported in get_stx_coordnum",
248 static void tpx_make_chain_identifiers(t_atoms
*atoms
, t_block
*mols
)
250 /* We always assign a new chain number, but save the chain id characters
251 * for larger molecules.
253 const int chainMinAtoms
= 15;
257 bool outOfIds
= false;
258 for (int m
= 0; m
< mols
->nr
; m
++)
260 int a0
= mols
->index
[m
];
261 int a1
= mols
->index
[m
+1];
263 if (a1
- a0
>= chainMinAtoms
&& !outOfIds
)
265 /* Set the chain id for the output */
267 /* Here we allow for the max possible 2*26+10=62 chain ids */
272 else if (chainid
== 'z')
276 else if (chainid
== '9')
289 for (int a
= a0
; a
< a1
; a
++)
291 atoms
->resinfo
[atoms
->atom
[a
].resind
].chainnum
= chainnum
;
292 atoms
->resinfo
[atoms
->atom
[a
].resind
].chainid
= c
;
297 /* Blank out the chain id if there was only one chain */
300 for (int r
= 0; r
< atoms
->nres
; r
++)
302 atoms
->resinfo
[r
].chainid
= ' ';
307 static void read_stx_conf(const char *infile
,
308 t_symtab
*symtab
, char ***name
, t_atoms
*atoms
,
309 rvec x
[], rvec
*v
, int *ePBC
, matrix box
)
314 char g96_line
[STRLEN
+1];
318 fprintf(stderr
, "Warning: Number of atoms in %s is 0\n", infile
);
320 else if (atoms
->atom
== NULL
)
322 gmx_mem("Uninitialized array atom");
330 ftp
= fn2ftp(infile
);
334 gmx_gro_read_conf(infile
, symtab
, name
, atoms
, x
, v
, box
);
338 fr
.natoms
= atoms
->nr
;
343 in
= gmx_fio_fopen(infile
, "r");
344 read_g96_conf(in
, infile
, &fr
, symtab
, g96_line
);
346 copy_mat(fr
.box
, box
);
347 *name
= put_symtab(symtab
, fr
.title
);
348 sfree(const_cast<char *>(fr
.title
));
353 gmx_pdb_read_conf(infile
, symtab
, name
, atoms
, x
, ePBC
, box
);
356 gmx_espresso_read_conf(infile
, symtab
, name
, atoms
, x
, v
, box
);
359 gmx_incons("Not supported in read_stx_conf");
363 static void done_gmx_groups_t(gmx_groups_t
*g
)
367 for (i
= 0; (i
< egcNR
); i
++)
369 if (NULL
!= g
->grps
[i
].nm_ind
)
371 sfree(g
->grps
[i
].nm_ind
);
372 g
->grps
[i
].nm_ind
= NULL
;
374 if (NULL
!= g
->grpnr
[i
])
380 /* The contents of this array is in symtab, don't free it here */
384 static void readConfAndAtoms(const char *infile
,
385 t_symtab
*symtab
, char ***name
, t_atoms
*atoms
,
387 rvec
**x
, rvec
**v
, matrix box
)
390 get_stx_coordnum(infile
, &natoms
);
392 init_t_atoms(atoms
, natoms
, (fn2ftp(infile
) == efPDB
));
394 bool xIsNull
= false;
405 read_stx_conf(infile
,
407 *x
, (v
== NULL
) ? NULL
: *v
, ePBC
, box
);
415 void readConfAndTopology(const char *infile
,
416 bool *haveTopology
, gmx_mtop_t
*mtop
,
418 rvec
**x
, rvec
**v
, matrix box
)
420 GMX_RELEASE_ASSERT(mtop
!= NULL
, "readConfAndTopology requires mtop!=NULL");
427 *haveTopology
= fn2bTPX(infile
);
431 read_tpxheader(infile
, &header
, TRUE
);
434 snew(*x
, header
.natoms
);
438 snew(*v
, header
.natoms
);
442 = read_tpx(infile
, NULL
, box
, &natoms
,
443 (x
== NULL
) ? NULL
: *x
, (v
== NULL
) ? NULL
: *v
, mtop
);
455 open_symtab(&symtab
);
457 readConfAndAtoms(infile
, &symtab
, &name
, &atoms
, ePBC
, x
, v
, box
);
460 convertAtomsToMtop(&symtab
, name
, &atoms
, mtop
);
464 gmx_bool
read_tps_conf(const char *infile
, t_topology
*top
, int *ePBC
,
465 rvec
**x
, rvec
**v
, matrix box
, gmx_bool requireMasses
)
470 // Note: We should have an initializer instead of relying on snew
472 readConfAndTopology(infile
, &haveTopology
, mtop
, ePBC
, x
, v
, box
);
474 *top
= gmx_mtop_t_to_t_topology(mtop
);
475 /* In this case we need to throw away the group data too */
476 done_gmx_groups_t(&mtop
->groups
);
479 tpx_make_chain_identifiers(&top
->atoms
, &top
->mols
);
481 if (requireMasses
&& !top
->atoms
.haveMass
)
483 atomsSetMassesBasedOnNames(&top
->atoms
, TRUE
);
485 if (!top
->atoms
.haveMass
)
487 gmx_fatal(FARGS
, "Masses were requested, but for some atom(s) masses could not be found in the database. Use a tpr file as input, if possible, or add these atoms to the mass database.");