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 main (int argc
,char *argv
[])
54 const char *desc
[] = {
55 "[TT]protonate[tt] reads (a) conformation(s) and adds all missing",
56 "hydrogens as defined in [TT]ffgmx2.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 pdb 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
;
80 int nidx
,natoms
,natoms_out
;
87 { efTPS
, NULL
, NULL
, ffREAD
},
88 { efTRX
, "-f", NULL
, ffOPTRD
},
89 { efNDX
, NULL
, NULL
, ffOPTRD
},
90 { efTRO
, "-o", "protonated", ffWRITE
}
92 #define NFILE asize(fnm)
94 CopyRight(stderr
,argv
[0]);
95 parse_common_args(&argc
,argv
,PCA_CAN_TIME
,
96 NFILE
,fnm
,0,NULL
,asize(desc
),desc
,0,NULL
,&oenv
);
98 infile
=opt2fn("-s",NFILE
,fnm
);
99 read_tps_conf(infile
,title
,&top
,&ePBC
,&x
,NULL
,box
,FALSE
);
101 printf("Select group to process:\n");
102 get_index(atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),1,&nidx
,&index
,&grpnm
);
103 bReadMultiple
= opt2bSet("-f",NFILE
,fnm
);
105 infile
= opt2fn("-f",NFILE
,fnm
);
106 if ( !read_first_frame(oenv
,&status
, infile
, &fr
, TRX_NEED_X
) )
107 gmx_fatal(FARGS
,"cannot read coordinate file %s",infile
);
110 clear_trxframe(&fr
,TRUE
);
111 fr
.natoms
= atoms
->nr
;
117 copy_mat(box
, fr
.box
);
123 gmx_fatal(FARGS
,"no atoms in coordinate file %s",infile
);
124 if ( natoms
> atoms
->nr
)
125 gmx_fatal(FARGS
,"topology with %d atoms does not match "
126 "coordinates with %d atoms",atoms
->nr
,natoms
);
127 for(i
=0; i
<nidx
; i
++)
128 if (index
[i
] > natoms
)
129 gmx_fatal(FARGS
,"An atom number in group %s is larger than the number of "
130 "atoms (%d) in the coordinate file %s",grpnm
,natoms
,infile
);
132 /* get indexed copy of atoms */
134 init_t_atoms(iatoms
,nidx
,FALSE
);
135 snew(iatoms
->atom
, iatoms
->nr
);
137 for(i
=0; i
<nidx
; i
++) {
138 iatoms
->atom
[i
] = atoms
->atom
[index
[i
]];
139 iatoms
->atomname
[i
] = atoms
->atomname
[index
[i
]];
140 if ( i
>0 && (atoms
->atom
[index
[i
]].resind
!=atoms
->atom
[index
[i
-1]].resind
) )
142 iatoms
->atom
[i
].resind
= resind
;
143 iatoms
->resinfo
[resind
] = atoms
->resinfo
[atoms
->atom
[index
[i
]].resind
];
144 iatoms
->nres
= max(iatoms
->nres
, iatoms
->atom
[i
].resind
+1);
147 init_t_protonate(&protdata
);
149 out
= open_trx(opt2fn("-o",NFILE
,fnm
),"w");
153 if (debug
) fprintf(debug
,"FRAME %d (%d %g)\n",frame
,fr
.step
,fr
.time
);
154 /* get indexed copy of x */
155 for(i
=0; i
<nidx
; i
++)
156 copy_rvec(fr
.x
[index
[i
]], ix
[i
]);
158 natoms_out
= protonate(&iatoms
, &ix
, &protdata
);
160 /* setup output frame */
162 frout
.natoms
= natoms_out
;
164 frout
.atoms
= iatoms
;
170 write_trxframe(out
,&frout
,NULL
);
172 } while ( bReadMultiple
&& read_next_frame(oenv
,status
, &fr
) );