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 main (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
);
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
);
111 infile
= opt2fn("-f",NFILE
,fnm
);
112 if ( !read_first_frame(oenv
,&status
, infile
, &fr
, TRX_NEED_X
) ) {
113 gmx_fatal(FARGS
,"cannot read coordinate file %s",infile
);
117 clear_trxframe(&fr
,TRUE
);
118 fr
.natoms
= atoms
->nr
;
124 copy_mat(box
, fr
.box
);
130 gmx_fatal(FARGS
,"no atoms in coordinate file %s",infile
);
133 if ( natoms
> atoms
->nr
) {
134 gmx_fatal(FARGS
,"topology with %d atoms does not match "
135 "coordinates with %d atoms",atoms
->nr
,natoms
);
138 for(i
=0; i
<nidx
; i
++) {
139 if (index
[i
] > natoms
) {
140 gmx_fatal(FARGS
,"An atom number in group %s is larger than the number of "
141 "atoms (%d) in the coordinate file %s",grpnm
,natoms
,infile
);
145 /* get indexed copy of atoms */
147 init_t_atoms(iatoms
,nidx
,FALSE
);
148 snew(iatoms
->atom
, iatoms
->nr
);
150 for(i
=0; i
<nidx
; i
++) {
151 iatoms
->atom
[i
] = atoms
->atom
[index
[i
]];
152 iatoms
->atomname
[i
] = atoms
->atomname
[index
[i
]];
153 if ( i
>0 && (atoms
->atom
[index
[i
]].resind
!=atoms
->atom
[index
[i
-1]].resind
) ) {
156 iatoms
->atom
[i
].resind
= resind
;
157 iatoms
->resinfo
[resind
] = atoms
->resinfo
[atoms
->atom
[index
[i
]].resind
];
158 /* allocate some space for the rtp name and copy from name */
159 snew(iatoms
->resinfo
[resind
].rtp
,1);
160 strcpy(iatoms
->resinfo
[resind
].rtp
, atoms
->resinfo
[resind
].name
);
162 iatoms
->nres
= max(iatoms
->nres
, iatoms
->atom
[i
].resind
+1);
165 init_t_protonate(&protdata
);
167 out
= open_trx(opt2fn("-o",NFILE
,fnm
),"w");
172 fprintf(debug
,"FRAME %d (%d %g)\n",frame
,fr
.step
,fr
.time
);
174 /* get indexed copy of x */
175 for(i
=0; i
<nidx
; i
++) {
176 copy_rvec(fr
.x
[index
[i
]], ix
[i
]);
179 natoms_out
= protonate(&iatoms
, &ix
, &protdata
);
181 /* setup output frame */
183 frout
.natoms
= natoms_out
;
185 frout
.atoms
= iatoms
;
191 write_trxframe(out
,&frout
,NULL
);
193 } while ( bReadMultiple
&& read_next_frame(oenv
,status
, &fr
) );