3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
51 #include "hackblock.h"
53 int cmain (int argc, char *argv[])
55 const char *desc[] = {
56 "[TT]g_protonate[tt] reads (a) conformation(s) and adds all missing",
57 "hydrogens as defined in [TT]gmx2.ff/aminoacids.hdb[tt]. If only [TT]-s[tt] is",
58 "specified, this conformation will be protonated, if also [TT]-f[tt]",
59 "is specified, the conformation(s) will be read from this file, ",
60 "which can be either a single conformation or a trajectory.",
62 "If a [TT].pdb[tt] file is supplied, residue names might not correspond to",
63 "to the GROMACS naming conventions, in which case these residues will",
64 "probably not be properly protonated.",
66 "If an index file is specified, please note that the atom numbers",
67 "should correspond to the [BB]protonated[bb] state."
75 t_atoms *atoms, *iatoms;
82 int nidx, natoms, natoms_out;
85 gmx_bool bReadMultiple;
88 const char *bugs[] = {
89 "For the moment, only .pdb files are accepted to the -s flag"
93 { efTPS, NULL, NULL, ffREAD },
94 { efTRX, "-f", NULL, ffOPTRD },
95 { efNDX, NULL, NULL, ffOPTRD },
96 { efTRO, "-o", "protonated", ffWRITE }
98 #define NFILE asize(fnm)
100 CopyRight(stderr, argv[0]);
101 parse_common_args(&argc, argv, PCA_CAN_TIME,
102 NFILE, fnm, 0, NULL, asize(desc), desc, asize(bugs), bugs, &oenv);
104 infile = opt2fn("-s", NFILE, fnm);
105 read_tps_conf(infile, title, &top, &ePBC, &x, NULL, box, FALSE);
106 atoms = &(top.atoms);
107 printf("Select group to process:\n");
108 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nidx, &index, &grpnm);
109 bReadMultiple = opt2bSet("-f", NFILE, fnm);
112 infile = opt2fn("-f", NFILE, fnm);
113 if (!read_first_frame(oenv, &status, infile, &fr, TRX_NEED_X ) )
115 gmx_fatal(FARGS, "cannot read coordinate file %s", infile);
121 clear_trxframe(&fr, TRUE);
122 fr.natoms = atoms->nr;
128 copy_mat(box, fr.box);
135 gmx_fatal(FARGS, "no atoms in coordinate file %s", infile);
138 if (natoms > atoms->nr)
140 gmx_fatal(FARGS, "topology with %d atoms does not match "
141 "coordinates with %d atoms", atoms->nr, natoms);
144 for (i = 0; i < nidx; i++)
146 if (index[i] > natoms)
148 gmx_fatal(FARGS, "An atom number in group %s is larger than the number of "
149 "atoms (%d) in the coordinate file %s", grpnm, natoms, infile);
153 /* get indexed copy of atoms */
155 init_t_atoms(iatoms, nidx, FALSE);
156 snew(iatoms->atom, iatoms->nr);
158 for (i = 0; i < nidx; i++)
160 iatoms->atom[i] = atoms->atom[index[i]];
161 iatoms->atomname[i] = atoms->atomname[index[i]];
162 if (i > 0 && (atoms->atom[index[i]].resind != atoms->atom[index[i-1]].resind) )
166 iatoms->atom[i].resind = resind;
167 iatoms->resinfo[resind] = atoms->resinfo[atoms->atom[index[i]].resind];
168 /* allocate some space for the rtp name and copy from name */
169 snew(iatoms->resinfo[resind].rtp, 1);
170 *iatoms->resinfo[resind].rtp = gmx_strdup(*atoms->resinfo[resind].name);
172 iatoms->nres = max(iatoms->nres, iatoms->atom[i].resind+1);
175 init_t_protonate(&protdata);
177 out = open_trx(opt2fn("-o", NFILE, fnm), "w");
184 fprintf(debug, "FRAME %d (%d %g)\n", frame, fr.step, fr.time);
186 /* get indexed copy of x */
187 for (i = 0; i < nidx; i++)
189 copy_rvec(fr.x[index[i]], ix[i]);
192 natoms_out = protonate(&iatoms, &ix, &protdata);
194 /* setup output frame */
196 frout.natoms = natoms_out;
198 frout.atoms = iatoms;
204 write_trxframe(out, &frout, NULL);
207 while (bReadMultiple && read_next_frame(oenv, status, &fr) );