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
50 #include "hackblock.h"
52 int gmx_protonate(int argc, char *argv[])
54 const char *desc[] = {
55 "[TT]g_protonate[tt] reads (a) conformation(s) and adds all missing",
56 "hydrogens as defined in [TT]gmx2.ff/aminoacids.hdb[tt]. If only [TT]-s[tt] is",
57 "specified, this conformation will be protonated, if also [TT]-f[tt]",
58 "is specified, the conformation(s) will be read from this file, ",
59 "which can be either a single conformation or a trajectory.",
61 "If a [TT].pdb[tt] file is supplied, residue names might not correspond to",
62 "to the GROMACS naming conventions, in which case these residues will",
63 "probably not be properly protonated.",
65 "If an index file is specified, please note that the atom numbers",
66 "should correspond to the [BB]protonated[bb] state."
74 t_atoms *atoms, *iatoms;
81 int nidx, natoms, natoms_out;
84 gmx_bool bReadMultiple;
87 const char *bugs[] = {
88 "For the moment, only .pdb files are accepted to the -s flag"
92 { efTPS, NULL, NULL, ffREAD },
93 { efTRX, "-f", NULL, ffOPTRD },
94 { efNDX, NULL, NULL, ffOPTRD },
95 { efTRO, "-o", "protonated", ffWRITE }
97 #define NFILE asize(fnm)
99 if (!parse_common_args(&argc, argv, PCA_CAN_TIME,
100 NFILE, fnm, 0, NULL, asize(desc), desc, asize(bugs), bugs, &oenv))
105 infile = opt2fn("-s", NFILE, fnm);
106 read_tps_conf(infile, title, &top, &ePBC, &x, NULL, box, FALSE);
107 atoms = &(top.atoms);
108 printf("Select group to process:\n");
109 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nidx, &index, &grpnm);
110 bReadMultiple = opt2bSet("-f", NFILE, fnm);
113 infile = opt2fn("-f", NFILE, fnm);
114 if (!read_first_frame(oenv, &status, infile, &fr, TRX_NEED_X ) )
116 gmx_fatal(FARGS, "cannot read coordinate file %s", infile);
122 clear_trxframe(&fr, TRUE);
123 fr.natoms = atoms->nr;
129 copy_mat(box, fr.box);
136 gmx_fatal(FARGS, "no atoms in coordinate file %s", infile);
139 if (natoms > atoms->nr)
141 gmx_fatal(FARGS, "topology with %d atoms does not match "
142 "coordinates with %d atoms", atoms->nr, natoms);
145 for (i = 0; i < nidx; i++)
147 if (index[i] > natoms)
149 gmx_fatal(FARGS, "An atom number in group %s is larger than the number of "
150 "atoms (%d) in the coordinate file %s", grpnm, natoms, infile);
154 /* get indexed copy of atoms */
156 init_t_atoms(iatoms, nidx, FALSE);
157 snew(iatoms->atom, iatoms->nr);
159 for (i = 0; i < nidx; i++)
161 iatoms->atom[i] = atoms->atom[index[i]];
162 iatoms->atomname[i] = atoms->atomname[index[i]];
163 if (i > 0 && (atoms->atom[index[i]].resind != atoms->atom[index[i-1]].resind) )
167 iatoms->atom[i].resind = resind;
168 iatoms->resinfo[resind] = atoms->resinfo[atoms->atom[index[i]].resind];
169 /* allocate some space for the rtp name and copy from name */
170 snew(iatoms->resinfo[resind].rtp, 1);
171 *iatoms->resinfo[resind].rtp = gmx_strdup(*atoms->resinfo[resind].name);
173 iatoms->nres = max(iatoms->nres, iatoms->atom[i].resind+1);
176 init_t_protonate(&protdata);
178 out = open_trx(opt2fn("-o", NFILE, fnm), "w");
185 fprintf(debug, "FRAME %d (%d %g)\n", frame, fr.step, fr.time);
187 /* get indexed copy of x */
188 for (i = 0; i < nidx; i++)
190 copy_rvec(fr.x[index[i]], ix[i]);
193 natoms_out = protonate(&iatoms, &ix, &protdata);
195 /* setup output frame */
197 frout.natoms = natoms_out;
199 frout.atoms = iatoms;
205 write_trxframe(out, &frout, NULL);
208 while (bReadMultiple && read_next_frame(oenv, status, &fr) );