Bug Summary

File:gromacs/gmxpreprocess/pdb2gmx.c
Location:line 1429, column 13
Description:Array access results in a null pointer dereference

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, 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.
10 *
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.
15 *
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.
20 *
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.
25 *
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.
33 *
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.
36 */
37#include "pdb2gmx.h"
38
39#ifdef HAVE_CONFIG_H1
40#include <config.h>
41#endif
42
43#include <ctype.h>
44#include <stdlib.h>
45#include <string.h>
46#include <time.h>
47
48#include "typedefs.h"
49#include "gromacs/fileio/gmxfio.h"
50#include "gromacs/utility/smalloc.h"
51#include "copyrite.h"
52#include "gromacs/utility/cstringutil.h"
53#include "gromacs/fileio/confio.h"
54#include "symtab.h"
55#include "gromacs/math/vec.h"
56#include "gromacs/commandline/pargs.h"
57#include "gromacs/utility/futil.h"
58#include "gromacs/utility/fatalerror.h"
59#include "gromacs/fileio/pdbio.h"
60#include "toputil.h"
61#include "h_db.h"
62#include "physics.h"
63#include "pgutil.h"
64#include "calch.h"
65#include "resall.h"
66#include "pdb2top.h"
67#include "ter_db.h"
68#include "gromacs/gmxlib/conformation-utilities.h"
69#include "genhydro.h"
70#include "readinp.h"
71#include "atomprop.h"
72#include "index.h"
73#include "fflibutil.h"
74#include "macros.h"
75
76#include "gromacs/fileio/strdb.h"
77
78#include "hizzie.h"
79#include "specbond.h"
80#include "xlate.h"
81
82typedef struct {
83 char gmx[6];
84 char main[6];
85 char nter[6];
86 char cter[6];
87 char bter[6];
88} rtprename_t;
89
90
91static const char *res2bb_notermini(const char *name,
92 int nrr, const rtprename_t *rr)
93{
94 /* NOTE: This function returns the main building block name,
95 * it does not take terminal renaming into account.
96 */
97 int i;
98
99 i = 0;
100 while (i < nrr && gmx_strcasecmp(name, rr[i].gmx) != 0)
101 {
102 i++;
103 }
104
105 return (i < nrr ? rr[i].main : name);
106}
107
108static const char *select_res(int nr, int resnr,
109 const char *name[], const char *expl[],
110 const char *title,
111 int nrr, const rtprename_t *rr)
112{
113 int sel = 0;
114
115 printf("Which %s type do you want for residue %d\n", title, resnr+1);
116 for (sel = 0; (sel < nr); sel++)
117 {
118 printf("%d. %s (%s)\n",
119 sel, expl[sel], res2bb_notermini(name[sel], nrr, rr));
120 }
121 printf("\nType a number:"); fflush(stdoutstdout);
122
123 if (scanf("%d", &sel) != 1)
124 {
125 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 125, "Answer me for res %s %d!", title, resnr+1);
126 }
127
128 return name[sel];
129}
130
131static const char *get_asptp(int resnr, int nrr, const rtprename_t *rr)
132{
133 enum {
134 easp, easpH, easpNR
135 };
136 const char *lh[easpNR] = { "ASP", "ASPH" };
137 const char *expl[easpNR] = {
138 "Not protonated (charge -1)",
139 "Protonated (charge 0)"
140 };
141
142 return select_res(easpNR, resnr, lh, expl, "ASPARTIC ACID", nrr, rr);
143}
144
145static const char *get_glutp(int resnr, int nrr, const rtprename_t *rr)
146{
147 enum {
148 eglu, egluH, egluNR
149 };
150 const char *lh[egluNR] = { "GLU", "GLUH" };
151 const char *expl[egluNR] = {
152 "Not protonated (charge -1)",
153 "Protonated (charge 0)"
154 };
155
156 return select_res(egluNR, resnr, lh, expl, "GLUTAMIC ACID", nrr, rr);
157}
158
159static const char *get_glntp(int resnr, int nrr, const rtprename_t *rr)
160{
161 enum {
162 egln, eglnH, eglnNR
163 };
164 const char *lh[eglnNR] = { "GLN", "QLN" };
165 const char *expl[eglnNR] = {
166 "Not protonated (charge 0)",
167 "Protonated (charge +1)"
168 };
169
170 return select_res(eglnNR, resnr, lh, expl, "GLUTAMINE", nrr, rr);
171}
172
173static const char *get_lystp(int resnr, int nrr, const rtprename_t *rr)
174{
175 enum {
176 elys, elysH, elysNR
177 };
178 const char *lh[elysNR] = { "LYSN", "LYS" };
179 const char *expl[elysNR] = {
180 "Not protonated (charge 0)",
181 "Protonated (charge +1)"
182 };
183
184 return select_res(elysNR, resnr, lh, expl, "LYSINE", nrr, rr);
185}
186
187static const char *get_argtp(int resnr, int nrr, const rtprename_t *rr)
188{
189 enum {
190 earg, eargH, eargNR
191 };
192 const char *lh[eargNR] = { "ARGN", "ARG" };
193 const char *expl[eargNR] = {
194 "Not protonated (charge 0)",
195 "Protonated (charge +1)"
196 };
197
198 return select_res(eargNR, resnr, lh, expl, "ARGININE", nrr, rr);
199}
200
201static const char *get_histp(int resnr, int nrr, const rtprename_t *rr)
202{
203 const char *expl[ehisNR] = {
204 "H on ND1 only",
205 "H on NE2 only",
206 "H on ND1 and NE2",
207 "Coupled to Heme"
208 };
209
210 return select_res(ehisNR, resnr, hh, expl, "HISTIDINE", nrr, rr);
211}
212
213static void read_rtprename(const char *fname, FILE *fp,
214 int *nrtprename, rtprename_t **rtprename)
215{
216 char line[STRLEN4096], buf[STRLEN4096];
217 int n;
218 rtprename_t *rr;
219 int ncol, nc;
220
221 n = *nrtprename;
222 rr = *rtprename;
223
224 ncol = 0;
225 while (get_a_line(fp, line, STRLEN4096))
226 {
227 srenew(rr, n+1)(rr) = save_realloc("rr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 227, (rr), (n+1), sizeof(*(rr)));
228 nc = sscanf(line, "%s %s %s %s %s %s",
229 rr[n].gmx, rr[n].main, rr[n].nter, rr[n].cter, rr[n].bter, buf);
230 if (ncol == 0)
231 {
232 if (nc != 2 && nc != 5)
233 {
234 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 234, "Residue renaming database '%s' has %d columns instead of %d, %d or %d", fname, ncol, 2, 5);
235 }
236 ncol = nc;
237 }
238 else if (nc != ncol)
239 {
240 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 240, "A line in residue renaming database '%s' has %d columns, while previous lines have %d columns", fname, nc, ncol);
241 }
242
243 if (nc == 2)
244 {
245 /* This file does not have special termini names, copy them from main */
246 strcpy(rr[n].nter, rr[n].main);
247 strcpy(rr[n].cter, rr[n].main);
248 strcpy(rr[n].bter, rr[n].main);
249 }
250
251 n++;
252 }
253
254 *nrtprename = n;
255 *rtprename = rr;
256}
257
258static char *search_resrename(int nrr, rtprename_t *rr,
259 const char *name,
260 gmx_bool bStart, gmx_bool bEnd,
261 gmx_bool bCompareFFRTPname)
262{
263 char *nn;
264 int i;
265
266 nn = NULL((void*)0);
267
268 i = 0;
269 while (i < nrr && ((!bCompareFFRTPname && strcmp(name, rr[i].gmx)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(name) && __builtin_constant_p (rr[i].gmx) &&
(__s1_len = strlen (name), __s2_len = strlen (rr[i].gmx), (!
((size_t)(const void *)((name) + 1) - (size_t)(const void *)(
name) == 1) || __s1_len >= 4) && (!((size_t)(const
void *)((rr[i].gmx) + 1) - (size_t)(const void *)(rr[i].gmx)
== 1) || __s2_len >= 4)) ? __builtin_strcmp (name, rr[i].
gmx) : (__builtin_constant_p (name) && ((size_t)(const
void *)((name) + 1) - (size_t)(const void *)(name) == 1) &&
(__s1_len = strlen (name), __s1_len < 4) ? (__builtin_constant_p
(rr[i].gmx) && ((size_t)(const void *)((rr[i].gmx) +
1) - (size_t)(const void *)(rr[i].gmx) == 1) ? __builtin_strcmp
(name, rr[i].gmx) : (__extension__ ({ const unsigned char *__s2
= (const unsigned char *) (const char *) (rr[i].gmx); int __result
= (((const unsigned char *) (const char *) (name))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (name))[1] - __s2
[1]); if (__s1_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (name))[2] - __s2
[2]); if (__s1_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) (name))[3] - __s2[3
]); } } __result; }))) : (__builtin_constant_p (rr[i].gmx) &&
((size_t)(const void *)((rr[i].gmx) + 1) - (size_t)(const void
*)(rr[i].gmx) == 1) && (__s2_len = strlen (rr[i].gmx
), __s2_len < 4) ? (__builtin_constant_p (name) &&
((size_t)(const void *)((name) + 1) - (size_t)(const void *)
(name) == 1) ? __builtin_strcmp (name, rr[i].gmx) : (- (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (name); int __result = (((const unsigned char *) (const
char *) (rr[i].gmx))[0] - __s2[0]); if (__s2_len > 0 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (rr[i].gmx))[1] - __s2[1]); if (__s2_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (rr[i].gmx))[2] - __s2[2]); if (__s2_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) (rr[i].gmx))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(name, rr[i].gmx)))); }) != 0) ||
270 ( bCompareFFRTPname && strcmp(name, rr[i].main)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(name) && __builtin_constant_p (rr[i].main) &&
(__s1_len = strlen (name), __s2_len = strlen (rr[i].main), (
!((size_t)(const void *)((name) + 1) - (size_t)(const void *)
(name) == 1) || __s1_len >= 4) && (!((size_t)(const
void *)((rr[i].main) + 1) - (size_t)(const void *)(rr[i].main
) == 1) || __s2_len >= 4)) ? __builtin_strcmp (name, rr[i]
.main) : (__builtin_constant_p (name) && ((size_t)(const
void *)((name) + 1) - (size_t)(const void *)(name) == 1) &&
(__s1_len = strlen (name), __s1_len < 4) ? (__builtin_constant_p
(rr[i].main) && ((size_t)(const void *)((rr[i].main)
+ 1) - (size_t)(const void *)(rr[i].main) == 1) ? __builtin_strcmp
(name, rr[i].main) : (__extension__ ({ const unsigned char *
__s2 = (const unsigned char *) (const char *) (rr[i].main); int
__result = (((const unsigned char *) (const char *) (name))[
0] - __s2[0]); if (__s1_len > 0 && __result == 0) {
__result = (((const unsigned char *) (const char *) (name))[
1] - __s2[1]); if (__s1_len > 1 && __result == 0) {
__result = (((const unsigned char *) (const char *) (name))[
2] - __s2[2]); if (__s1_len > 2 && __result == 0) __result
= (((const unsigned char *) (const char *) (name))[3] - __s2
[3]); } } __result; }))) : (__builtin_constant_p (rr[i].main)
&& ((size_t)(const void *)((rr[i].main) + 1) - (size_t
)(const void *)(rr[i].main) == 1) && (__s2_len = strlen
(rr[i].main), __s2_len < 4) ? (__builtin_constant_p (name
) && ((size_t)(const void *)((name) + 1) - (size_t)(const
void *)(name) == 1) ? __builtin_strcmp (name, rr[i].main) : (
- (__extension__ ({ const unsigned char *__s2 = (const unsigned
char *) (const char *) (name); int __result = (((const unsigned
char *) (const char *) (rr[i].main))[0] - __s2[0]); if (__s2_len
> 0 && __result == 0) { __result = (((const unsigned
char *) (const char *) (rr[i].main))[1] - __s2[1]); if (__s2_len
> 1 && __result == 0) { __result = (((const unsigned
char *) (const char *) (rr[i].main))[2] - __s2[2]); if (__s2_len
> 2 && __result == 0) __result = (((const unsigned
char *) (const char *) (rr[i].main))[3] - __s2[3]); } } __result
; })))) : __builtin_strcmp (name, rr[i].main)))); }) != 0)))
271 {
272 i++;
273 }
274
275 /* If found in the database, rename this residue's rtp building block,
276 * otherwise keep the old name.
277 */
278 if (i < nrr)
279 {
280 if (bStart && bEnd)
281 {
282 nn = rr[i].bter;
283 }
284 else if (bStart)
285 {
286 nn = rr[i].nter;
287 }
288 else if (bEnd)
289 {
290 nn = rr[i].cter;
291 }
292 else
293 {
294 nn = rr[i].main;
295 }
296
297 if (nn[0] == '-')
298 {
299 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 299, "In the chosen force field there is no residue type for '%s'%s", name, bStart ? ( bEnd ? " as a standalone (starting & ending) residue" : " as a starting terminus") : (bEnd ? " as an ending terminus" : ""));
300 }
301 }
302
303 return nn;
304}
305
306
307static void rename_resrtp(t_atoms *pdba, int nterpairs, int *r_start, int *r_end,
308 int nrr, rtprename_t *rr, t_symtab *symtab,
309 gmx_bool bVerbose)
310{
311 int r, i, j;
312 gmx_bool bStart, bEnd;
313 char *nn;
314 gmx_bool bFFRTPTERRNM;
315
316 bFFRTPTERRNM = (getenv("GMX_NO_FFRTP_TER_RENAME") == NULL((void*)0));
317
318 for (r = 0; r < pdba->nres; r++)
319 {
320 bStart = FALSE0;
321 bEnd = FALSE0;
322 for (j = 0; j < nterpairs; j++)
323 {
324 if (r == r_start[j])
325 {
326 bStart = TRUE1;
327 }
328 }
329 for (j = 0; j < nterpairs; j++)
330 {
331 if (r == r_end[j])
332 {
333 bEnd = TRUE1;
334 }
335 }
336
337 nn = search_resrename(nrr, rr, *pdba->resinfo[r].rtp, bStart, bEnd, FALSE0);
338
339 if (bFFRTPTERRNM && nn == NULL((void*)0) && (bStart || bEnd))
340 {
341 /* This is a terminal residue, but the residue name,
342 * currently stored in .rtp, is not a standard residue name,
343 * but probably a force field specific rtp name.
344 * Check if we need to rename it because it is terminal.
345 */
346 nn = search_resrename(nrr, rr,
347 *pdba->resinfo[r].rtp, bStart, bEnd, TRUE1);
348 }
349
350 if (nn != NULL((void*)0) && strcmp(*pdba->resinfo[r].rtp, nn)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(*pdba->resinfo[r].rtp) && __builtin_constant_p (
nn) && (__s1_len = strlen (*pdba->resinfo[r].rtp),
__s2_len = strlen (nn), (!((size_t)(const void *)((*pdba->
resinfo[r].rtp) + 1) - (size_t)(const void *)(*pdba->resinfo
[r].rtp) == 1) || __s1_len >= 4) && (!((size_t)(const
void *)((nn) + 1) - (size_t)(const void *)(nn) == 1) || __s2_len
>= 4)) ? __builtin_strcmp (*pdba->resinfo[r].rtp, nn) :
(__builtin_constant_p (*pdba->resinfo[r].rtp) && (
(size_t)(const void *)((*pdba->resinfo[r].rtp) + 1) - (size_t
)(const void *)(*pdba->resinfo[r].rtp) == 1) && (__s1_len
= strlen (*pdba->resinfo[r].rtp), __s1_len < 4) ? (__builtin_constant_p
(nn) && ((size_t)(const void *)((nn) + 1) - (size_t)
(const void *)(nn) == 1) ? __builtin_strcmp (*pdba->resinfo
[r].rtp, nn) : (__extension__ ({ const unsigned char *__s2 = (
const unsigned char *) (const char *) (nn); int __result = ((
(const unsigned char *) (const char *) (*pdba->resinfo[r].
rtp))[0] - __s2[0]); if (__s1_len > 0 && __result ==
0) { __result = (((const unsigned char *) (const char *) (*pdba
->resinfo[r].rtp))[1] - __s2[1]); if (__s1_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (*pdba->resinfo[r].rtp))[2] - __s2[2]); if (__s1_len
> 2 && __result == 0) __result = (((const unsigned
char *) (const char *) (*pdba->resinfo[r].rtp))[3] - __s2
[3]); } } __result; }))) : (__builtin_constant_p (nn) &&
((size_t)(const void *)((nn) + 1) - (size_t)(const void *)(nn
) == 1) && (__s2_len = strlen (nn), __s2_len < 4) ?
(__builtin_constant_p (*pdba->resinfo[r].rtp) && (
(size_t)(const void *)((*pdba->resinfo[r].rtp) + 1) - (size_t
)(const void *)(*pdba->resinfo[r].rtp) == 1) ? __builtin_strcmp
(*pdba->resinfo[r].rtp, nn) : (- (__extension__ ({ const unsigned
char *__s2 = (const unsigned char *) (const char *) (*pdba->
resinfo[r].rtp); int __result = (((const unsigned char *) (const
char *) (nn))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
nn))[1] - __s2[1]); if (__s2_len > 1 && __result ==
0) { __result = (((const unsigned char *) (const char *) (nn
))[2] - __s2[2]); if (__s2_len > 2 && __result == 0
) __result = (((const unsigned char *) (const char *) (nn))[3
] - __s2[3]); } } __result; })))) : __builtin_strcmp (*pdba->
resinfo[r].rtp, nn)))); }) != 0)
351 {
352 if (bVerbose)
353 {
354 printf("Changing rtp entry of residue %d %s to '%s'\n",
355 pdba->resinfo[r].nr, *pdba->resinfo[r].name, nn);
356 }
357 pdba->resinfo[r].rtp = put_symtab(symtab, nn);
358 }
359 }
360}
361
362static void pdbres_to_gmxrtp(t_atoms *pdba)
363{
364 int i;
365
366 for (i = 0; (i < pdba->nres); i++)
367 {
368 if (pdba->resinfo[i].rtp == NULL((void*)0))
369 {
370 pdba->resinfo[i].rtp = pdba->resinfo[i].name;
371 }
372 }
373}
374
375static void rename_pdbres(t_atoms *pdba, const char *oldnm, const char *newnm,
376 gmx_bool bFullCompare, t_symtab *symtab)
377{
378 char *resnm;
379 int i;
380
381 for (i = 0; (i < pdba->nres); i++)
382 {
383 resnm = *pdba->resinfo[i].name;
384 if ((bFullCompare && (gmx_strcasecmp(resnm, oldnm) == 0)) ||
385 (!bFullCompare && strstr(resnm, oldnm) != NULL((void*)0)))
386 {
387 /* Rename the residue name (not the rtp name) */
388 pdba->resinfo[i].name = put_symtab(symtab, newnm);
389 }
390 }
391}
392
393static void rename_bb(t_atoms *pdba, const char *oldnm, const char *newnm,
394 gmx_bool bFullCompare, t_symtab *symtab)
395{
396 char *bbnm;
397 int i;
398
399 for (i = 0; (i < pdba->nres); i++)
400 {
401 /* We have not set the rtp name yes, use the residue name */
402 bbnm = *pdba->resinfo[i].name;
403 if ((bFullCompare && (gmx_strcasecmp(bbnm, oldnm) == 0)) ||
404 (!bFullCompare && strstr(bbnm, oldnm) != NULL((void*)0)))
405 {
406 /* Change the rtp builing block name */
407 pdba->resinfo[i].rtp = put_symtab(symtab, newnm);
408 }
409 }
410}
411
412static void rename_bbint(t_atoms *pdba, const char *oldnm,
413 const char *gettp(int, int, const rtprename_t *),
414 gmx_bool bFullCompare,
415 t_symtab *symtab,
416 int nrr, const rtprename_t *rr)
417{
418 int i;
419 const char *ptr;
420 char *bbnm;
421
422 for (i = 0; i < pdba->nres; i++)
423 {
424 /* We have not set the rtp name yes, use the residue name */
425 bbnm = *pdba->resinfo[i].name;
426 if ((bFullCompare && (strcmp(bbnm, oldnm)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(bbnm) && __builtin_constant_p (oldnm) && (__s1_len
= strlen (bbnm), __s2_len = strlen (oldnm), (!((size_t)(const
void *)((bbnm) + 1) - (size_t)(const void *)(bbnm) == 1) || __s1_len
>= 4) && (!((size_t)(const void *)((oldnm) + 1) -
(size_t)(const void *)(oldnm) == 1) || __s2_len >= 4)) ? __builtin_strcmp
(bbnm, oldnm) : (__builtin_constant_p (bbnm) && ((size_t
)(const void *)((bbnm) + 1) - (size_t)(const void *)(bbnm) ==
1) && (__s1_len = strlen (bbnm), __s1_len < 4) ? (
__builtin_constant_p (oldnm) && ((size_t)(const void *
)((oldnm) + 1) - (size_t)(const void *)(oldnm) == 1) ? __builtin_strcmp
(bbnm, oldnm) : (__extension__ ({ const unsigned char *__s2 =
(const unsigned char *) (const char *) (oldnm); int __result
= (((const unsigned char *) (const char *) (bbnm))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (bbnm))[1] - __s2
[1]); if (__s1_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (bbnm))[2] - __s2
[2]); if (__s1_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) (bbnm))[3] - __s2[3
]); } } __result; }))) : (__builtin_constant_p (oldnm) &&
((size_t)(const void *)((oldnm) + 1) - (size_t)(const void *
)(oldnm) == 1) && (__s2_len = strlen (oldnm), __s2_len
< 4) ? (__builtin_constant_p (bbnm) && ((size_t)(
const void *)((bbnm) + 1) - (size_t)(const void *)(bbnm) == 1
) ? __builtin_strcmp (bbnm, oldnm) : (- (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
(bbnm); int __result = (((const unsigned char *) (const char
*) (oldnm))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
oldnm))[1] - __s2[1]); if (__s2_len > 1 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
oldnm))[2] - __s2[2]); if (__s2_len > 2 && __result
== 0) __result = (((const unsigned char *) (const char *) (oldnm
))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp (bbnm
, oldnm)))); }) == 0)) ||
427 (!bFullCompare && strstr(bbnm, oldnm) != NULL((void*)0)))
428 {
429 ptr = gettp(i, nrr, rr);
430 pdba->resinfo[i].rtp = put_symtab(symtab, ptr);
431 }
432 }
433}
434
435static void check_occupancy(t_atoms *atoms, const char *filename, gmx_bool bVerbose)
436{
437 int i, ftp;
438 int nzero = 0;
439 int nnotone = 0;
440
441 ftp = fn2ftp(filename);
442 if (!atoms->pdbinfo || ((ftp != efPDB) && (ftp != efBRK) && (ftp != efENT)))
443 {
444 fprintf(stderrstderr, "No occupancies in %s\n", filename);
445 }
446 else
447 {
448 for (i = 0; (i < atoms->nr); i++)
449 {
450 if (atoms->pdbinfo[i].occup != 1)
451 {
452 if (bVerbose)
453 {
454 fprintf(stderrstderr, "Occupancy for atom %s%d-%s is %f rather than 1\n",
455 *atoms->resinfo[atoms->atom[i].resind].name,
456 atoms->resinfo[ atoms->atom[i].resind].nr,
457 *atoms->atomname[i],
458 atoms->pdbinfo[i].occup);
459 }
460 if (atoms->pdbinfo[i].occup == 0)
461 {
462 nzero++;
463 }
464 else
465 {
466 nnotone++;
467 }
468 }
469 }
470 if (nzero == atoms->nr)
471 {
472 fprintf(stderrstderr, "All occupancy fields zero. This is probably not an X-Ray structure\n");
473 }
474 else if ((nzero > 0) || (nnotone > 0))
475 {
476 fprintf(stderrstderr,
477 "\n"
478 "WARNING: there were %d atoms with zero occupancy and %d atoms with\n"
479 " occupancy unequal to one (out of %d atoms). Check your pdb file.\n"
480 "\n",
481 nzero, nnotone, atoms->nr);
482 }
483 else
484 {
485 fprintf(stderrstderr, "All occupancies are one\n");
486 }
487 }
488}
489
490void write_posres(char *fn, t_atoms *pdba, real fc)
491{
492 FILE *fp;
493 int i;
494
495 fp = gmx_fio_fopen(fn, "w");
496 fprintf(fp,
497 "; In this topology include file, you will find position restraint\n"
498 "; entries for all the heavy atoms in your original pdb file.\n"
499 "; This means that all the protons which were added by pdb2gmx are\n"
500 "; not restrained.\n"
501 "\n"
502 "[ position_restraints ]\n"
503 "; %4s%6s%8s%8s%8s\n", "atom", "type", "fx", "fy", "fz"
504 );
505 for (i = 0; (i < pdba->nr); i++)
506 {
507 if (!is_hydrogen(*pdba->atomname[i]) && !is_dummymass(*pdba->atomname[i]))
508 {
509 fprintf(fp, "%6d%6d %g %g %g\n", i+1, 1, fc, fc, fc);
510 }
511 }
512 gmx_fio_fclose(fp);
513}
514
515static int read_pdball(const char *inf, const char *outf, char *title,
516 t_atoms *atoms, rvec **x,
517 int *ePBC, matrix box, gmx_bool bRemoveH,
518 t_symtab *symtab, gmx_residuetype_t rt, const char *watres,
519 gmx_atomprop_t aps, gmx_bool bVerbose)
520/* Read a pdb file. (containing proteins) */
521{
522 int natom, new_natom, i;
523
524 /* READ IT */
525 printf("Reading %s...\n", inf);
526 get_stx_coordnum(inf, &natom);
527 init_t_atoms(atoms, natom, TRUE1);
528 snew(*x, natom)(*x) = save_calloc("*x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 528, (natom), sizeof(*(*x)));
529 read_stx_conf(inf, title, atoms, *x, NULL((void*)0), ePBC, box);
530 if (fn2ftp(inf) == efPDB)
531 {
532 get_pdb_atomnumber(atoms, aps);
533 }
534 if (bRemoveH)
535 {
536 new_natom = 0;
537 for (i = 0; i < atoms->nr; i++)
538 {
539 if (!is_hydrogen(*atoms->atomname[i]))
540 {
541 atoms->atom[new_natom] = atoms->atom[i];
542 atoms->atomname[new_natom] = atoms->atomname[i];
543 atoms->pdbinfo[new_natom] = atoms->pdbinfo[i];
544 copy_rvec((*x)[i], (*x)[new_natom]);
545 new_natom++;
546 }
547 }
548 atoms->nr = new_natom;
549 natom = new_natom;
550 }
551
552 printf("Read");
553 if (title && title[0])
554 {
555 printf(" '%s',", title);
556 }
557 printf(" %d atoms\n", natom);
558
559 /* Rename residues */
560 rename_pdbres(atoms, "HOH", watres, FALSE0, symtab);
561 rename_pdbres(atoms, "SOL", watres, FALSE0, symtab);
562 rename_pdbres(atoms, "WAT", watres, FALSE0, symtab);
563
564 rename_atoms("xlateat.dat", NULL((void*)0),
565 atoms, symtab, NULL((void*)0), TRUE1, rt, TRUE1, bVerbose);
566
567 if (natom == 0)
568 {
569 return 0;
570 }
571
572 if (outf)
573 {
574 write_sto_conf(outf, title, atoms, *x, NULL((void*)0), *ePBC, box);
575 }
576
577 return natom;
578}
579
580void process_chain(t_atoms *pdba, rvec *x,
581 gmx_bool bTrpU, gmx_bool bPheU, gmx_bool bTyrU,
582 gmx_bool bLysMan, gmx_bool bAspMan, gmx_bool bGluMan,
583 gmx_bool bHisMan, gmx_bool bArgMan, gmx_bool bGlnMan,
584 real angle, real distance, t_symtab *symtab,
585 int nrr, const rtprename_t *rr)
586{
587 /* Rename aromatics, lys, asp and histidine */
588 if (bTyrU)
589 {
590 rename_bb(pdba, "TYR", "TYRU", FALSE0, symtab);
591 }
592 if (bTrpU)
593 {
594 rename_bb(pdba, "TRP", "TRPU", FALSE0, symtab);
595 }
596 if (bPheU)
597 {
598 rename_bb(pdba, "PHE", "PHEU", FALSE0, symtab);
599 }
600 if (bLysMan)
601 {
602 rename_bbint(pdba, "LYS", get_lystp, FALSE0, symtab, nrr, rr);
603 }
604 if (bArgMan)
605 {
606 rename_bbint(pdba, "ARG", get_argtp, FALSE0, symtab, nrr, rr);
607 }
608 if (bGlnMan)
609 {
610 rename_bbint(pdba, "GLN", get_glntp, FALSE0, symtab, nrr, rr);
611 }
612 if (bAspMan)
613 {
614 rename_bbint(pdba, "ASP", get_asptp, FALSE0, symtab, nrr, rr);
615 }
616 else
617 {
618 rename_bb(pdba, "ASPH", "ASP", FALSE0, symtab);
619 }
620 if (bGluMan)
621 {
622 rename_bbint(pdba, "GLU", get_glutp, FALSE0, symtab, nrr, rr);
623 }
624 else
625 {
626 rename_bb(pdba, "GLUH", "GLU", FALSE0, symtab);
627 }
628
629 if (!bHisMan)
630 {
631 set_histp(pdba, x, angle, distance);
632 }
633 else
634 {
635 rename_bbint(pdba, "HIS", get_histp, TRUE1, symtab, nrr, rr);
636 }
637
638 /* Initialize the rtp builing block names with the residue names
639 * for the residues that have not been processed above.
640 */
641 pdbres_to_gmxrtp(pdba);
642
643 /* Now we have all rtp names set.
644 * The rtp names will conform to Gromacs naming,
645 * unless the input pdb file contained one or more force field specific
646 * rtp names as residue names.
647 */
648}
649
650/* struct for sorting the atoms from the pdb file */
651typedef struct {
652 int resnr; /* residue number */
653 int j; /* database order index */
654 int index; /* original atom number */
655 char anm1; /* second letter of atom name */
656 char altloc; /* alternate location indicator */
657} t_pdbindex;
658
659int pdbicomp(const void *a, const void *b)
660{
661 t_pdbindex *pa, *pb;
662 int d;
663
664 pa = (t_pdbindex *)a;
665 pb = (t_pdbindex *)b;
666
667 d = (pa->resnr - pb->resnr);
668 if (d == 0)
669 {
670 d = (pa->j - pb->j);
671 if (d == 0)
672 {
673 d = (pa->anm1 - pb->anm1);
674 if (d == 0)
675 {
676 d = (pa->altloc - pb->altloc);
677 }
678 }
679 }
680
681 return d;
682}
683
684static void sort_pdbatoms(t_restp restp[],
685 int natoms, t_atoms **pdbaptr, rvec **x,
686 t_blocka *block, char ***gnames)
687{
688 t_atoms *pdba, *pdbnew;
689 rvec **xnew;
690 int i, j;
691 t_restp *rptr;
692 t_hackblock *hbr;
693 t_pdbindex *pdbi;
694 atom_id *a;
695 char *atomnm;
696
697 pdba = *pdbaptr;
698 natoms = pdba->nr;
699 pdbnew = NULL((void*)0);
700 snew(xnew, 1)(xnew) = save_calloc("xnew", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 700, (1), sizeof(*(xnew)));
701 snew(pdbi, natoms)(pdbi) = save_calloc("pdbi", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 701, (natoms), sizeof(*(pdbi)));
702
703 for (i = 0; i < natoms; i++)
704 {
705 atomnm = *pdba->atomname[i];
706 rptr = &restp[pdba->atom[i].resind];
707 for (j = 0; (j < rptr->natom); j++)
708 {
709 if (gmx_strcasecmp(atomnm, *(rptr->atomname[j])) == 0)
710 {
711 break;
712 }
713 }
714 if (j == rptr->natom)
715 {
716 char buf[STRLEN4096];
717
718 sprintf(buf,
719 "Atom %s in residue %s %d was not found in rtp entry %s with %d atoms\n"
720 "while sorting atoms.\n%s", atomnm,
721 *pdba->resinfo[pdba->atom[i].resind].name,
722 pdba->resinfo[pdba->atom[i].resind].nr,
723 rptr->resname,
724 rptr->natom,
725 is_hydrogen(atomnm) ?
726 "\nFor a hydrogen, this can be a different protonation state, or it\n"
727 "might have had a different number in the PDB file and was rebuilt\n"
728 "(it might for instance have been H3, and we only expected H1 & H2).\n"
729 "Note that hydrogens might have been added to the entry for the N-terminus.\n"
730 "Remove this hydrogen or choose a different protonation state to solve it.\n"
731 "Option -ignh will ignore all hydrogens in the input." : ".");
732 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 732, buf);
733 }
734 /* make shadow array to be sorted into indexgroup */
735 pdbi[i].resnr = pdba->atom[i].resind;
736 pdbi[i].j = j;
737 pdbi[i].index = i;
738 pdbi[i].anm1 = atomnm[1];
739 pdbi[i].altloc = pdba->pdbinfo[i].altloc;
740 }
741 qsort(pdbi, natoms, (size_t)sizeof(pdbi[0]), pdbicomp);
742
743 /* pdba is sorted in pdbnew using the pdbi index */
744 snew(a, natoms)(a) = save_calloc("a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 744, (natoms), sizeof(*(a)));
745 snew(pdbnew, 1)(pdbnew) = save_calloc("pdbnew", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 745, (1), sizeof(*(pdbnew)));
746 init_t_atoms(pdbnew, natoms, TRUE1);
747 snew(*xnew, natoms)(*xnew) = save_calloc("*xnew", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 747, (natoms), sizeof(*(*xnew)));
748 pdbnew->nr = pdba->nr;
749 pdbnew->nres = pdba->nres;
750 sfree(pdbnew->resinfo)save_free("pdbnew->resinfo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 750, (pdbnew->resinfo));
751 pdbnew->resinfo = pdba->resinfo;
752 for (i = 0; i < natoms; i++)
753 {
754 pdbnew->atom[i] = pdba->atom[pdbi[i].index];
755 pdbnew->atomname[i] = pdba->atomname[pdbi[i].index];
756 pdbnew->pdbinfo[i] = pdba->pdbinfo[pdbi[i].index];
757 copy_rvec((*x)[pdbi[i].index], (*xnew)[i]);
758 /* make indexgroup in block */
759 a[i] = pdbi[i].index;
760 }
761 /* clean up */
762 sfree(pdba->atomname)save_free("pdba->atomname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 762, (pdba->atomname));
763 sfree(pdba->atom)save_free("pdba->atom", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 763, (pdba->atom));
764 sfree(pdba->pdbinfo)save_free("pdba->pdbinfo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 764, (pdba->pdbinfo));
765 sfree(pdba)save_free("pdba", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 765, (pdba));
766 sfree(*x)save_free("*x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 766, (*x));
767 /* copy the sorted pdbnew back to pdba */
768 *pdbaptr = pdbnew;
769 *x = *xnew;
770 add_grp(block, gnames, natoms, a, "prot_sort");
771 sfree(xnew)save_free("xnew", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 771, (xnew));
772 sfree(a)save_free("a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 772, (a));
773 sfree(pdbi)save_free("pdbi", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 773, (pdbi));
774}
775
776static int remove_duplicate_atoms(t_atoms *pdba, rvec x[], gmx_bool bVerbose)
777{
778 int i, j, oldnatoms, ndel;
779 t_resinfo *ri;
780
781 printf("Checking for duplicate atoms....\n");
782 oldnatoms = pdba->nr;
783 ndel = 0;
784 /* NOTE: pdba->nr is modified inside the loop */
785 for (i = 1; (i < pdba->nr); i++)
786 {
787 /* compare 'i' and 'i-1', throw away 'i' if they are identical
788 this is a 'while' because multiple alternate locations can be present */
789 while ( (i < pdba->nr) &&
790 (pdba->atom[i-1].resind == pdba->atom[i].resind) &&
791 (strcmp(*pdba->atomname[i-1], *pdba->atomname[i])__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(*pdba->atomname[i-1]) && __builtin_constant_p (*
pdba->atomname[i]) && (__s1_len = strlen (*pdba->
atomname[i-1]), __s2_len = strlen (*pdba->atomname[i]), (!
((size_t)(const void *)((*pdba->atomname[i-1]) + 1) - (size_t
)(const void *)(*pdba->atomname[i-1]) == 1) || __s1_len >=
4) && (!((size_t)(const void *)((*pdba->atomname[
i]) + 1) - (size_t)(const void *)(*pdba->atomname[i]) == 1
) || __s2_len >= 4)) ? __builtin_strcmp (*pdba->atomname
[i-1], *pdba->atomname[i]) : (__builtin_constant_p (*pdba->
atomname[i-1]) && ((size_t)(const void *)((*pdba->
atomname[i-1]) + 1) - (size_t)(const void *)(*pdba->atomname
[i-1]) == 1) && (__s1_len = strlen (*pdba->atomname
[i-1]), __s1_len < 4) ? (__builtin_constant_p (*pdba->atomname
[i]) && ((size_t)(const void *)((*pdba->atomname[i
]) + 1) - (size_t)(const void *)(*pdba->atomname[i]) == 1)
? __builtin_strcmp (*pdba->atomname[i-1], *pdba->atomname
[i]) : (__extension__ ({ const unsigned char *__s2 = (const unsigned
char *) (const char *) (*pdba->atomname[i]); int __result
= (((const unsigned char *) (const char *) (*pdba->atomname
[i-1]))[0] - __s2[0]); if (__s1_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
*pdba->atomname[i-1]))[1] - __s2[1]); if (__s1_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (*pdba->atomname[i-1]))[2] - __s2[2]); if (__s1_len
> 2 && __result == 0) __result = (((const unsigned
char *) (const char *) (*pdba->atomname[i-1]))[3] - __s2[
3]); } } __result; }))) : (__builtin_constant_p (*pdba->atomname
[i]) && ((size_t)(const void *)((*pdba->atomname[i
]) + 1) - (size_t)(const void *)(*pdba->atomname[i]) == 1)
&& (__s2_len = strlen (*pdba->atomname[i]), __s2_len
< 4) ? (__builtin_constant_p (*pdba->atomname[i-1]) &&
((size_t)(const void *)((*pdba->atomname[i-1]) + 1) - (size_t
)(const void *)(*pdba->atomname[i-1]) == 1) ? __builtin_strcmp
(*pdba->atomname[i-1], *pdba->atomname[i]) : (- (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (*pdba->atomname[i-1]); int __result = (((const unsigned
char *) (const char *) (*pdba->atomname[i]))[0] - __s2[0]
); if (__s2_len > 0 && __result == 0) { __result =
(((const unsigned char *) (const char *) (*pdba->atomname
[i]))[1] - __s2[1]); if (__s2_len > 1 && __result ==
0) { __result = (((const unsigned char *) (const char *) (*pdba
->atomname[i]))[2] - __s2[2]); if (__s2_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) (*pdba->atomname[i]))[3] - __s2[3]); } } __result; }))
)) : __builtin_strcmp (*pdba->atomname[i-1], *pdba->atomname
[i])))); }) == 0) )
792 {
793 ndel++;
794 if (bVerbose)
795 {
796 ri = &pdba->resinfo[pdba->atom[i].resind];
797 printf("deleting duplicate atom %4s %s%4d%c",
798 *pdba->atomname[i], *ri->name, ri->nr, ri->ic);
799 if (ri->chainid && (ri->chainid != ' '))
800 {
801 printf(" ch %c", ri->chainid);
802 }
803 if (pdba->pdbinfo)
804 {
805 if (pdba->pdbinfo[i].atomnr)
806 {
807 printf(" pdb nr %4d", pdba->pdbinfo[i].atomnr);
808 }
809 if (pdba->pdbinfo[i].altloc && (pdba->pdbinfo[i].altloc != ' '))
810 {
811 printf(" altloc %c", pdba->pdbinfo[i].altloc);
812 }
813 }
814 printf("\n");
815 }
816 pdba->nr--;
817 /* We can not free, since it might be in the symtab */
818 /* sfree(pdba->atomname[i]); */
819 for (j = i; j < pdba->nr; j++)
820 {
821 pdba->atom[j] = pdba->atom[j+1];
822 pdba->atomname[j] = pdba->atomname[j+1];
823 pdba->pdbinfo[j] = pdba->pdbinfo[j+1];
824 copy_rvec(x[j+1], x[j]);
825 }
826 srenew(pdba->atom, pdba->nr)(pdba->atom) = save_realloc("pdba->atom", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 826, (pdba->atom), (pdba->nr), sizeof(*(pdba->atom
)));
827 /* srenew(pdba->atomname, pdba->nr); */
828 srenew(pdba->pdbinfo, pdba->nr)(pdba->pdbinfo) = save_realloc("pdba->pdbinfo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 828, (pdba->pdbinfo), (pdba->nr), sizeof(*(pdba->pdbinfo
)));
829 }
830 }
831 if (pdba->nr != oldnatoms)
832 {
833 printf("Now there are %d atoms. Deleted %d duplicates.\n", pdba->nr, ndel);
834 }
835
836 return pdba->nr;
837}
838
839void find_nc_ter(t_atoms *pdba, int r0, int r1, int *r_start, int *r_end, gmx_residuetype_t rt)
840{
841 int i;
842 const char *p_startrestype;
843 const char *p_restype;
844 int nstartwarn, nendwarn;
845
846 *r_start = -1;
847 *r_end = -1;
848
849 nstartwarn = 0;
850 nendwarn = 0;
851
852 /* Find the starting terminus (typially N or 5') */
853 for (i = r0; i < r1 && *r_start == -1; i++)
854 {
855 gmx_residuetype_get_type(rt, *pdba->resinfo[i].name, &p_startrestype);
856 if (!gmx_strcasecmp(p_startrestype, "Protein") || !gmx_strcasecmp(p_startrestype, "DNA") || !gmx_strcasecmp(p_startrestype, "RNA") )
857 {
858 printf("Identified residue %s%d as a starting terminus.\n", *pdba->resinfo[i].name, pdba->resinfo[i].nr);
859 *r_start = i;
860 }
861 else
862 {
863 if (nstartwarn < 5)
864 {
865 printf("Warning: Starting residue %s%d in chain not identified as Protein/RNA/DNA.\n", *pdba->resinfo[i].name, pdba->resinfo[i].nr);
866 }
867 if (nstartwarn == 5)
868 {
869 printf("More than 5 unidentified residues at start of chain - disabling further warnings.\n");
870 }
871 nstartwarn++;
872 }
873 }
874
875 if (*r_start >= 0)
876 {
877 /* Go through the rest of the residues, check that they are the same class, and identify the ending terminus. */
878 for (i = *r_start; i < r1; i++)
879 {
880 gmx_residuetype_get_type(rt, *pdba->resinfo[i].name, &p_restype);
881 if (!gmx_strcasecmp(p_restype, p_startrestype) && nendwarn == 0)
882 {
883 *r_end = i;
884 }
885 else
886 {
887 if (nendwarn < 5)
888 {
889 printf("Warning: Residue %s%d in chain has different type (%s) from starting residue %s%d (%s).\n",
890 *pdba->resinfo[i].name, pdba->resinfo[i].nr, p_restype,
891 *pdba->resinfo[*r_start].name, pdba->resinfo[*r_start].nr, p_startrestype);
892 }
893 if (nendwarn == 5)
894 {
895 printf("More than 5 unidentified residues at end of chain - disabling further warnings.\n");
896 }
897 nendwarn++;
898 }
899 }
900 }
901
902 if (*r_end >= 0)
903 {
904 printf("Identified residue %s%d as a ending terminus.\n", *pdba->resinfo[*r_end].name, pdba->resinfo[*r_end].nr);
905 }
906}
907
908
909static void
910modify_chain_numbers(t_atoms * pdba,
911 const char * chainsep)
912{
913 int i;
914 char old_prev_chainid;
915 char old_this_chainid;
916 int old_prev_chainnum;
917 int old_this_chainnum;
918 t_resinfo *ri;
919 char select[STRLEN4096];
920 int new_chainnum;
921 int this_atomnum;
922 int prev_atomnum;
923 const char * prev_atomname;
924 const char * this_atomname;
925 const char * prev_resname;
926 const char * this_resname;
927 int prev_resnum;
928 int this_resnum;
929 char prev_chainid;
930 char this_chainid;
931 int prev_chainnumber;
932 int this_chainnumber;
933
934 enum
935 {
936 SPLIT_ID_OR_TER,
937 SPLIT_ID_AND_TER,
938 SPLIT_ID_ONLY,
939 SPLIT_TER_ONLY,
940 SPLIT_INTERACTIVE
941 }
942 splitting;
943
944 splitting = SPLIT_TER_ONLY; /* keep compiler happy */
945
946 /* Be a bit flexible to catch typos */
947 if (!strncmp(chainsep, "id_o", 4)(__extension__ (__builtin_constant_p (4) && ((__builtin_constant_p
(chainsep) && strlen (chainsep) < ((size_t) (4)))
|| (__builtin_constant_p ("id_o") && strlen ("id_o")
< ((size_t) (4)))) ? __extension__ ({ size_t __s1_len, __s2_len
; (__builtin_constant_p (chainsep) && __builtin_constant_p
("id_o") && (__s1_len = strlen (chainsep), __s2_len =
strlen ("id_o"), (!((size_t)(const void *)((chainsep) + 1) -
(size_t)(const void *)(chainsep) == 1) || __s1_len >= 4) &&
(!((size_t)(const void *)(("id_o") + 1) - (size_t)(const void
*)("id_o") == 1) || __s2_len >= 4)) ? __builtin_strcmp (chainsep
, "id_o") : (__builtin_constant_p (chainsep) && ((size_t
)(const void *)((chainsep) + 1) - (size_t)(const void *)(chainsep
) == 1) && (__s1_len = strlen (chainsep), __s1_len <
4) ? (__builtin_constant_p ("id_o") && ((size_t)(const
void *)(("id_o") + 1) - (size_t)(const void *)("id_o") == 1)
? __builtin_strcmp (chainsep, "id_o") : (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
("id_o"); int __result = (((const unsigned char *) (const char
*) (chainsep))[0] - __s2[0]); if (__s1_len > 0 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (chainsep))[1] - __s2[1]); if (__s1_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (chainsep))[2] - __s2[2]); if (__s1_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) (chainsep))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
("id_o") && ((size_t)(const void *)(("id_o") + 1) - (
size_t)(const void *)("id_o") == 1) && (__s2_len = strlen
("id_o"), __s2_len < 4) ? (__builtin_constant_p (chainsep
) && ((size_t)(const void *)((chainsep) + 1) - (size_t
)(const void *)(chainsep) == 1) ? __builtin_strcmp (chainsep,
"id_o") : (- (__extension__ ({ const unsigned char *__s2 = (
const unsigned char *) (const char *) (chainsep); int __result
= (((const unsigned char *) (const char *) ("id_o"))[0] - __s2
[0]); if (__s2_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) ("id_o"))[1] - __s2
[1]); if (__s2_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) ("id_o"))[2] - __s2
[2]); if (__s2_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) ("id_o"))[3] - __s2
[3]); } } __result; })))) : __builtin_strcmp (chainsep, "id_o"
)))); }) : strncmp (chainsep, "id_o", 4))))
948 {
949 /* For later interactive splitting we tentatively assign new chain numbers at either changing id or ter records */
950 splitting = SPLIT_ID_OR_TER;
951 printf("Splitting chemical chains based on TER records or chain id changing.\n");
952 }
953 else if (!strncmp(chainsep, "int", 3)(__extension__ (__builtin_constant_p (3) && ((__builtin_constant_p
(chainsep) && strlen (chainsep) < ((size_t) (3)))
|| (__builtin_constant_p ("int") && strlen ("int") <
((size_t) (3)))) ? __extension__ ({ size_t __s1_len, __s2_len
; (__builtin_constant_p (chainsep) && __builtin_constant_p
("int") && (__s1_len = strlen (chainsep), __s2_len =
strlen ("int"), (!((size_t)(const void *)((chainsep) + 1) - (
size_t)(const void *)(chainsep) == 1) || __s1_len >= 4) &&
(!((size_t)(const void *)(("int") + 1) - (size_t)(const void
*)("int") == 1) || __s2_len >= 4)) ? __builtin_strcmp (chainsep
, "int") : (__builtin_constant_p (chainsep) && ((size_t
)(const void *)((chainsep) + 1) - (size_t)(const void *)(chainsep
) == 1) && (__s1_len = strlen (chainsep), __s1_len <
4) ? (__builtin_constant_p ("int") && ((size_t)(const
void *)(("int") + 1) - (size_t)(const void *)("int") == 1) ?
__builtin_strcmp (chainsep, "int") : (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
("int"); int __result = (((const unsigned char *) (const char
*) (chainsep))[0] - __s2[0]); if (__s1_len > 0 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (chainsep))[1] - __s2[1]); if (__s1_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (chainsep))[2] - __s2[2]); if (__s1_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) (chainsep))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
("int") && ((size_t)(const void *)(("int") + 1) - (size_t
)(const void *)("int") == 1) && (__s2_len = strlen ("int"
), __s2_len < 4) ? (__builtin_constant_p (chainsep) &&
((size_t)(const void *)((chainsep) + 1) - (size_t)(const void
*)(chainsep) == 1) ? __builtin_strcmp (chainsep, "int") : (-
(__extension__ ({ const unsigned char *__s2 = (const unsigned
char *) (const char *) (chainsep); int __result = (((const unsigned
char *) (const char *) ("int"))[0] - __s2[0]); if (__s2_len >
0 && __result == 0) { __result = (((const unsigned char
*) (const char *) ("int"))[1] - __s2[1]); if (__s2_len > 1
&& __result == 0) { __result = (((const unsigned char
*) (const char *) ("int"))[2] - __s2[2]); if (__s2_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) ("int"))[3] - __s2[3]); } } __result; })))) :
__builtin_strcmp (chainsep, "int")))); }) : strncmp (chainsep
, "int", 3))))
954 {
955 /* For later interactive splitting we tentatively assign new chain numbers at either changing id or ter records */
956 splitting = SPLIT_INTERACTIVE;
957 printf("Splitting chemical chains interactively.\n");
958 }
959 else if (!strncmp(chainsep, "id_a", 4)(__extension__ (__builtin_constant_p (4) && ((__builtin_constant_p
(chainsep) && strlen (chainsep) < ((size_t) (4)))
|| (__builtin_constant_p ("id_a") && strlen ("id_a")
< ((size_t) (4)))) ? __extension__ ({ size_t __s1_len, __s2_len
; (__builtin_constant_p (chainsep) && __builtin_constant_p
("id_a") && (__s1_len = strlen (chainsep), __s2_len =
strlen ("id_a"), (!((size_t)(const void *)((chainsep) + 1) -
(size_t)(const void *)(chainsep) == 1) || __s1_len >= 4) &&
(!((size_t)(const void *)(("id_a") + 1) - (size_t)(const void
*)("id_a") == 1) || __s2_len >= 4)) ? __builtin_strcmp (chainsep
, "id_a") : (__builtin_constant_p (chainsep) && ((size_t
)(const void *)((chainsep) + 1) - (size_t)(const void *)(chainsep
) == 1) && (__s1_len = strlen (chainsep), __s1_len <
4) ? (__builtin_constant_p ("id_a") && ((size_t)(const
void *)(("id_a") + 1) - (size_t)(const void *)("id_a") == 1)
? __builtin_strcmp (chainsep, "id_a") : (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
("id_a"); int __result = (((const unsigned char *) (const char
*) (chainsep))[0] - __s2[0]); if (__s1_len > 0 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (chainsep))[1] - __s2[1]); if (__s1_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (chainsep))[2] - __s2[2]); if (__s1_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) (chainsep))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
("id_a") && ((size_t)(const void *)(("id_a") + 1) - (
size_t)(const void *)("id_a") == 1) && (__s2_len = strlen
("id_a"), __s2_len < 4) ? (__builtin_constant_p (chainsep
) && ((size_t)(const void *)((chainsep) + 1) - (size_t
)(const void *)(chainsep) == 1) ? __builtin_strcmp (chainsep,
"id_a") : (- (__extension__ ({ const unsigned char *__s2 = (
const unsigned char *) (const char *) (chainsep); int __result
= (((const unsigned char *) (const char *) ("id_a"))[0] - __s2
[0]); if (__s2_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) ("id_a"))[1] - __s2
[1]); if (__s2_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) ("id_a"))[2] - __s2
[2]); if (__s2_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) ("id_a"))[3] - __s2
[3]); } } __result; })))) : __builtin_strcmp (chainsep, "id_a"
)))); }) : strncmp (chainsep, "id_a", 4))))
960 {
961 splitting = SPLIT_ID_AND_TER;
962 printf("Splitting chemical chains based on TER records and chain id changing.\n");
963 }
964 else if (strlen(chainsep) == 2 && !strncmp(chainsep, "id", 4)(__extension__ (__builtin_constant_p (4) && ((__builtin_constant_p
(chainsep) && strlen (chainsep) < ((size_t) (4)))
|| (__builtin_constant_p ("id") && strlen ("id") <
((size_t) (4)))) ? __extension__ ({ size_t __s1_len, __s2_len
; (__builtin_constant_p (chainsep) && __builtin_constant_p
("id") && (__s1_len = strlen (chainsep), __s2_len = strlen
("id"), (!((size_t)(const void *)((chainsep) + 1) - (size_t)
(const void *)(chainsep) == 1) || __s1_len >= 4) &&
(!((size_t)(const void *)(("id") + 1) - (size_t)(const void *
)("id") == 1) || __s2_len >= 4)) ? __builtin_strcmp (chainsep
, "id") : (__builtin_constant_p (chainsep) && ((size_t
)(const void *)((chainsep) + 1) - (size_t)(const void *)(chainsep
) == 1) && (__s1_len = strlen (chainsep), __s1_len <
4) ? (__builtin_constant_p ("id") && ((size_t)(const
void *)(("id") + 1) - (size_t)(const void *)("id") == 1) ? __builtin_strcmp
(chainsep, "id") : (__extension__ ({ const unsigned char *__s2
= (const unsigned char *) (const char *) ("id"); int __result
= (((const unsigned char *) (const char *) (chainsep))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) (chainsep))[1] - __s2
[1]); if (__s1_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) (chainsep))[2] - __s2
[2]); if (__s1_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) (chainsep))[3] - __s2
[3]); } } __result; }))) : (__builtin_constant_p ("id") &&
((size_t)(const void *)(("id") + 1) - (size_t)(const void *)
("id") == 1) && (__s2_len = strlen ("id"), __s2_len <
4) ? (__builtin_constant_p (chainsep) && ((size_t)(const
void *)((chainsep) + 1) - (size_t)(const void *)(chainsep) ==
1) ? __builtin_strcmp (chainsep, "id") : (- (__extension__ (
{ const unsigned char *__s2 = (const unsigned char *) (const char
*) (chainsep); int __result = (((const unsigned char *) (const
char *) ("id"))[0] - __s2[0]); if (__s2_len > 0 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) ("id"))[1] - __s2[1]); if (__s2_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) ("id"))[2] - __s2[2]); if (__s2_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) ("id"))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(chainsep, "id")))); }) : strncmp (chainsep, "id", 4))))
965 {
966 splitting = SPLIT_ID_ONLY;
967 printf("Splitting chemical chains based on changing chain id only (ignoring TER records).\n");
968 }
969 else if (chainsep[0] == 't')
970 {
971 splitting = SPLIT_TER_ONLY;
972 printf("Splitting chemical chains based on TER records only (ignoring chain id).\n");
973 }
974 else
975 {
976 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 976, "Unidentified setting for chain separation: %s\n", chainsep);
977 }
978
979 /* The default chain enumeration is based on TER records only, which is reflected in chainnum below */
980
981 old_prev_chainid = '?';
982 old_prev_chainnum = -1;
983 new_chainnum = -1;
984
985 this_atomname = NULL((void*)0);
986 this_atomnum = -1;
987 this_resname = NULL((void*)0);
988 this_resnum = -1;
989 this_chainid = '?';
990 this_chainnumber = -1;
991
992 for (i = 0; i < pdba->nres; i++)
993 {
994 ri = &pdba->resinfo[i];
995 old_this_chainid = ri->chainid;
996 old_this_chainnum = ri->chainnum;
997
998 prev_atomname = this_atomname;
999 prev_atomnum = this_atomnum;
1000 prev_resname = this_resname;
1001 prev_resnum = this_resnum;
1002 prev_chainid = this_chainid;
1003 prev_chainnumber = this_chainnumber;
1004
1005 this_atomname = *(pdba->atomname[i]);
1006 this_atomnum = (pdba->pdbinfo != NULL((void*)0)) ? pdba->pdbinfo[i].atomnr : i+1;
1007 this_resname = *ri->name;
1008 this_resnum = ri->nr;
1009 this_chainid = ri->chainid;
1010 this_chainnumber = ri->chainnum;
1011
1012 switch (splitting)
1013 {
1014 case SPLIT_ID_OR_TER:
1015 if (old_this_chainid != old_prev_chainid || old_this_chainnum != old_prev_chainnum)
1016 {
1017 new_chainnum++;
1018 }
1019 break;
1020
1021 case SPLIT_ID_AND_TER:
1022 if (old_this_chainid != old_prev_chainid && old_this_chainnum != old_prev_chainnum)
1023 {
1024 new_chainnum++;
1025 }
1026 break;
1027
1028 case SPLIT_ID_ONLY:
1029 if (old_this_chainid != old_prev_chainid)
1030 {
1031 new_chainnum++;
1032 }
1033 break;
1034
1035 case SPLIT_TER_ONLY:
1036 if (old_this_chainnum != old_prev_chainnum)
1037 {
1038 new_chainnum++;
1039 }
1040 break;
1041 case SPLIT_INTERACTIVE:
1042 if (old_this_chainid != old_prev_chainid || old_this_chainnum != old_prev_chainnum)
1043 {
1044 if (i > 0)
1045 {
1046 printf("Split the chain (and introduce termini) between residue %s%d (chain id '%c', atom %d %s)\n"
1047 "and residue %s%d (chain id '%c', atom %d %s) ? [n/y]\n",
1048 prev_resname, prev_resnum, prev_chainid, prev_atomnum, prev_atomname,
1049 this_resname, this_resnum, this_chainid, this_atomnum, this_atomname);
1050
1051 if (NULL((void*)0) == fgets(select, STRLEN4096-1, stdinstdin))
1052 {
1053 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1053, "Error reading from stdin");
1054 }
1055 }
1056 if (i == 0 || select[0] == 'y')
1057 {
1058 new_chainnum++;
1059 }
1060 }
1061 break;
1062 default:
1063 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1063, "Internal inconsistency - this shouldn't happen...");
1064 break;
1065 }
1066 old_prev_chainid = old_this_chainid;
1067 old_prev_chainnum = old_this_chainnum;
1068
1069 ri->chainnum = new_chainnum;
1070 }
1071}
1072
1073
1074typedef struct {
1075 char chainid;
1076 char chainnum;
1077 int start;
1078 int natom;
1079 gmx_bool bAllWat;
1080 int nterpairs;
1081 int *chainstart;
1082} t_pdbchain;
1083
1084typedef struct {
1085 char chainid;
1086 int chainnum;
1087 gmx_bool bAllWat;
1088 int nterpairs;
1089 int *chainstart;
1090 t_hackblock **ntdb;
1091 t_hackblock **ctdb;
1092 int *r_start;
1093 int *r_end;
1094 t_atoms *pdba;
1095 rvec *x;
1096} t_chain;
1097
1098int gmx_pdb2gmx(int argc, char *argv[])
1099{
1100 const char *desc[] = {
1101 "[THISMODULE] reads a [TT].pdb[tt] (or [TT].gro[tt]) file, reads",
1102 "some database files, adds hydrogens to the molecules and generates",
1103 "coordinates in GROMACS (GROMOS), or optionally [TT].pdb[tt], format",
1104 "and a topology in GROMACS format.",
1105 "These files can subsequently be processed to generate a run input file.",
1106 "[PAR]",
1107 "[THISMODULE] will search for force fields by looking for",
1108 "a [TT]forcefield.itp[tt] file in subdirectories [TT]<forcefield>.ff[tt]",
1109 "of the current working directory and of the GROMACS library directory",
1110 "as inferred from the path of the binary or the [TT]GMXLIB[tt] environment",
1111 "variable.",
1112 "By default the forcefield selection is interactive,",
1113 "but you can use the [TT]-ff[tt] option to specify one of the short names",
1114 "in the list on the command line instead. In that case [THISMODULE] just looks",
1115 "for the corresponding [TT]<forcefield>.ff[tt] directory.",
1116 "[PAR]",
1117 "After choosing a force field, all files will be read only from",
1118 "the corresponding force field directory.",
1119 "If you want to modify or add a residue types, you can copy the force",
1120 "field directory from the GROMACS library directory to your current",
1121 "working directory. If you want to add new protein residue types,",
1122 "you will need to modify [TT]residuetypes.dat[tt] in the library directory",
1123 "or copy the whole library directory to a local directory and set",
1124 "the environment variable [TT]GMXLIB[tt] to the name of that directory.",
1125 "Check Chapter 5 of the manual for more information about file formats.",
1126 "[PAR]",
1127
1128 "Note that a [TT].pdb[tt] file is nothing more than a file format, and it",
1129 "need not necessarily contain a protein structure. Every kind of",
1130 "molecule for which there is support in the database can be converted.",
1131 "If there is no support in the database, you can add it yourself.[PAR]",
1132
1133 "The program has limited intelligence, it reads a number of database",
1134 "files, that allow it to make special bonds (Cys-Cys, Heme-His, etc.),",
1135 "if necessary this can be done manually. The program can prompt the",
1136 "user to select which kind of LYS, ASP, GLU, CYS or HIS residue is",
1137 "desired. For Lys the choice is between neutral (two protons on NZ) or",
1138 "protonated (three protons, default), for Asp and Glu unprotonated",
1139 "(default) or protonated, for His the proton can be either on ND1,",
1140 "on NE2 or on both. By default these selections are done automatically.",
1141 "For His, this is based on an optimal hydrogen bonding",
1142 "conformation. Hydrogen bonds are defined based on a simple geometric",
1143 "criterion, specified by the maximum hydrogen-donor-acceptor angle",
1144 "and donor-acceptor distance, which are set by [TT]-angle[tt] and",
1145 "[TT]-dist[tt] respectively.[PAR]",
1146
1147 "The protonation state of N- and C-termini can be chosen interactively",
1148 "with the [TT]-ter[tt] flag. Default termini are ionized (NH3+ and COO-),",
1149 "respectively. Some force fields support zwitterionic forms for chains of",
1150 "one residue, but for polypeptides these options should NOT be selected.",
1151 "The AMBER force fields have unique forms for the terminal residues,",
1152 "and these are incompatible with the [TT]-ter[tt] mechanism. You need",
1153 "to prefix your N- or C-terminal residue names with \"N\" or \"C\"",
1154 "respectively to use these forms, making sure you preserve the format",
1155 "of the coordinate file. Alternatively, use named terminating residues",
1156 "(e.g. ACE, NME).[PAR]",
1157
1158 "The separation of chains is not entirely trivial since the markup",
1159 "in user-generated PDB files frequently varies and sometimes it",
1160 "is desirable to merge entries across a TER record, for instance",
1161 "if you want a disulfide bridge or distance restraints between",
1162 "two protein chains or if you have a HEME group bound to a protein.",
1163 "In such cases multiple chains should be contained in a single",
1164 "[TT]moleculetype[tt] definition.",
1165 "To handle this, [THISMODULE] uses two separate options.",
1166 "First, [TT]-chainsep[tt] allows you to choose when a new chemical chain should",
1167 "start, and termini added when applicable. This can be done based on the",
1168 "existence of TER records, when the chain id changes, or combinations of either",
1169 "or both of these. You can also do the selection fully interactively.",
1170 "In addition, there is a [TT]-merge[tt] option that controls how multiple chains",
1171 "are merged into one moleculetype, after adding all the chemical termini (or not).",
1172 "This can be turned off (no merging), all non-water chains can be merged into a",
1173 "single molecule, or the selection can be done interactively.[PAR]",
1174
1175 "[THISMODULE] will also check the occupancy field of the [TT].pdb[tt] file.",
1176 "If any of the occupancies are not one, indicating that the atom is",
1177 "not resolved well in the structure, a warning message is issued.",
1178 "When a [TT].pdb[tt] file does not originate from an X-ray structure determination",
1179 "all occupancy fields may be zero. Either way, it is up to the user",
1180 "to verify the correctness of the input data (read the article!).[PAR]",
1181
1182 "During processing the atoms will be reordered according to GROMACS",
1183 "conventions. With [TT]-n[tt] an index file can be generated that",
1184 "contains one group reordered in the same way. This allows you to",
1185 "convert a GROMOS trajectory and coordinate file to GROMOS. There is",
1186 "one limitation: reordering is done after the hydrogens are stripped",
1187 "from the input and before new hydrogens are added. This means that",
1188 "you should not use [TT]-ignh[tt].[PAR]",
1189
1190 "The [TT].gro[tt] and [TT].g96[tt] file formats do not support chain",
1191 "identifiers. Therefore it is useful to enter a [TT].pdb[tt] file name at",
1192 "the [TT]-o[tt] option when you want to convert a multi-chain [TT].pdb[tt] file.",
1193 "[PAR]",
1194
1195 "The option [TT]-vsite[tt] removes hydrogen and fast improper dihedral",
1196 "motions. Angular and out-of-plane motions can be removed by changing",
1197 "hydrogens into virtual sites and fixing angles, which fixes their",
1198 "position relative to neighboring atoms. Additionally, all atoms in the",
1199 "aromatic rings of the standard amino acids (i.e. PHE, TRP, TYR and HIS)",
1200 "can be converted into virtual sites, eliminating the fast improper dihedral",
1201 "fluctuations in these rings. [BB]Note[bb] that in this case all other hydrogen",
1202 "atoms are also converted to virtual sites. The mass of all atoms that are",
1203 "converted into virtual sites, is added to the heavy atoms.[PAR]",
1204 "Also slowing down of dihedral motion can be done with [TT]-heavyh[tt]",
1205 "done by increasing the hydrogen-mass by a factor of 4. This is also",
1206 "done for water hydrogens to slow down the rotational motion of water.",
1207 "The increase in mass of the hydrogens is subtracted from the bonded",
1208 "(heavy) atom so that the total mass of the system remains the same."
1209 };
1210
1211
1212 FILE *fp, *top_file, *top_file2, *itp_file = NULL((void*)0);
1213 int natom, nres;
1214 t_atoms pdba_all, *pdba;
1215 t_atoms *atoms;
1216 t_resinfo *ri;
1217 t_blocka *block;
1218 int chain, nch, maxch, nwaterchain;
1219 t_pdbchain *pdb_ch;
1220 t_chain *chains, *cc;
1221 char select[STRLEN4096];
1222 int nincl, nmol;
1223 char **incls;
1224 t_mols *mols;
1225 char **gnames;
1226 int ePBC;
1227 matrix box;
1228 rvec box_space;
1229 int i, j, k, l, nrtp;
1230 int *swap_index, si;
1231 t_restp *restp;
1232 t_hackblock *ah;
1233 t_symtab symtab;
1234 gpp_atomtype_t atype;
1235 gmx_residuetype_t rt;
1236 const char *top_fn;
1237 char fn[256], itp_fn[STRLEN4096], posre_fn[STRLEN4096], buf_fn[STRLEN4096];
1238 char molname[STRLEN4096], title[STRLEN4096], quote[STRLEN4096];
1239 char *c, forcefield[STRLEN4096], ffdir[STRLEN4096];
1240 char ffname[STRLEN4096], suffix[STRLEN4096], buf[STRLEN4096];
1241 char *watermodel;
1242 const char *watres;
1243 int nrtpf;
1244 char **rtpf;
1245 char rtp[STRLEN4096];
1246 int nrrn;
1247 char **rrn;
1248 int nrtprename, naa;
1249 rtprename_t *rtprename = NULL((void*)0);
1250 int nah, nNtdb, nCtdb, ntdblist;
1251 t_hackblock *ntdb, *ctdb, **tdblist;
1252 int nssbonds;
1253 t_ssbond *ssbonds;
1254 rvec *pdbx, *x;
1255 gmx_bool bVsites = FALSE0, bWat, bPrevWat = FALSE0, bITP, bVsiteAromatics = FALSE0, bCheckMerge;
1256 real mHmult = 0;
1257 t_hackblock *hb_chain;
1258 t_restp *restp_chain;
1259 output_env_t oenv;
1260 const char *p_restype;
1261 int rc;
1262 int this_atomnum;
1263 int prev_atomnum;
1264 const char * prev_atomname;
1265 const char * this_atomname;
1266 const char * prev_resname;
1267 const char * this_resname;
1268 int prev_resnum;
1269 int this_resnum;
1270 char prev_chainid;
1271 char this_chainid;
1272 int prev_chainnumber;
1273 int this_chainnumber;
1274 int nid_used;
1275 int this_chainstart;
1276 int prev_chainstart;
1277 gmx_bool bMerged;
1278 int nchainmerges;
1279
1280 gmx_atomprop_t aps;
1281
1282 t_filenm fnm[] = {
1283 { efSTX, "-f", "eiwit.pdb", ffREAD1<<1 },
1284 { efSTO, "-o", "conf", ffWRITE1<<2 },
1285 { efTOP, NULL((void*)0), NULL((void*)0), ffWRITE1<<2 },
1286 { efITP, "-i", "posre", ffWRITE1<<2 },
1287 { efNDX, "-n", "clean", ffOPTWR(1<<2| 1<<3) },
1288 { efSTO, "-q", "clean.pdb", ffOPTWR(1<<2| 1<<3) }
1289 };
1290#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
1291
1292
1293 /* Command line arguments must be static */
1294 static gmx_bool bNewRTP = FALSE0;
1295 static gmx_bool bInter = FALSE0, bCysMan = FALSE0;
1296 static gmx_bool bLysMan = FALSE0, bAspMan = FALSE0, bGluMan = FALSE0, bHisMan = FALSE0;
1297 static gmx_bool bGlnMan = FALSE0, bArgMan = FALSE0;
1298 static gmx_bool bTerMan = FALSE0, bUnA = FALSE0, bHeavyH;
1299 static gmx_bool bSort = TRUE1, bAllowMissing = FALSE0, bRemoveH = FALSE0;
1300 static gmx_bool bDeuterate = FALSE0, bVerbose = FALSE0, bChargeGroups = TRUE1, bCmap = TRUE1;
1301 static gmx_bool bRenumRes = FALSE0, bRTPresname = FALSE0;
1302 static real angle = 135.0, distance = 0.3, posre_fc = 1000;
1303 static real long_bond_dist = 0.25, short_bond_dist = 0.05;
1304 static const char *vsitestr[] = { NULL((void*)0), "none", "hydrogens", "aromatics", NULL((void*)0) };
1305 static const char *watstr[] = { NULL((void*)0), "select", "none", "spc", "spce", "tip3p", "tip4p", "tip5p", NULL((void*)0) };
1306 static const char *chainsep[] = { NULL((void*)0), "id_or_ter", "id_and_ter", "ter", "id", "interactive", NULL((void*)0) };
1307 static const char *merge[] = {NULL((void*)0), "no", "all", "interactive", NULL((void*)0) };
1308 static const char *ff = "select";
1309
1310 t_pargs pa[] = {
1311 { "-newrtp", FALSE0, etBOOL, {&bNewRTP},
1312 "HIDDENWrite the residue database in new format to [TT]new.rtp[tt]"},
1313 { "-lb", FALSE0, etREAL, {&long_bond_dist},
1314 "HIDDENLong bond warning distance" },
1315 { "-sb", FALSE0, etREAL, {&short_bond_dist},
1316 "HIDDENShort bond warning distance" },
1317 { "-chainsep", FALSE0, etENUM, {chainsep},
1318 "Condition in PDB files when a new chain should be started (adding termini)" },
1319 { "-merge", FALSE0, etENUM, {&merge},
1320 "Merge multiple chains into a single [moleculetype]" },
1321 { "-ff", FALSE0, etSTR, {&ff},
1322 "Force field, interactive by default. Use [TT]-h[tt] for information." },
1323 { "-water", FALSE0, etENUM, {watstr},
1324 "Water model to use" },
1325 { "-inter", FALSE0, etBOOL, {&bInter},
1326 "Set the next 8 options to interactive"},
1327 { "-ss", FALSE0, etBOOL, {&bCysMan},
1328 "Interactive SS bridge selection" },
1329 { "-ter", FALSE0, etBOOL, {&bTerMan},
1330 "Interactive termini selection, instead of charged (default)" },
1331 { "-lys", FALSE0, etBOOL, {&bLysMan},
1332 "Interactive lysine selection, instead of charged" },
1333 { "-arg", FALSE0, etBOOL, {&bArgMan},
1334 "Interactive arginine selection, instead of charged" },
1335 { "-asp", FALSE0, etBOOL, {&bAspMan},
1336 "Interactive aspartic acid selection, instead of charged" },
1337 { "-glu", FALSE0, etBOOL, {&bGluMan},
1338 "Interactive glutamic acid selection, instead of charged" },
1339 { "-gln", FALSE0, etBOOL, {&bGlnMan},
1340 "Interactive glutamine selection, instead of neutral" },
1341 { "-his", FALSE0, etBOOL, {&bHisMan},
1342 "Interactive histidine selection, instead of checking H-bonds" },
1343 { "-angle", FALSE0, etREAL, {&angle},
1344 "Minimum hydrogen-donor-acceptor angle for a H-bond (degrees)" },
1345 { "-dist", FALSE0, etREAL, {&distance},
1346 "Maximum donor-acceptor distance for a H-bond (nm)" },
1347 { "-una", FALSE0, etBOOL, {&bUnA},
1348 "Select aromatic rings with united CH atoms on phenylalanine, "
1349 "tryptophane and tyrosine" },
1350 { "-sort", FALSE0, etBOOL, {&bSort},
1351 "HIDDENSort the residues according to database, turning this off is dangerous as charge groups might be broken in parts" },
1352 { "-ignh", FALSE0, etBOOL, {&bRemoveH},
1353 "Ignore hydrogen atoms that are in the coordinate file" },
1354 { "-missing", FALSE0, etBOOL, {&bAllowMissing},
1355 "Continue when atoms are missing, dangerous" },
1356 { "-v", FALSE0, etBOOL, {&bVerbose},
1357 "Be slightly more verbose in messages" },
1358 { "-posrefc", FALSE0, etREAL, {&posre_fc},
1359 "Force constant for position restraints" },
1360 { "-vsite", FALSE0, etENUM, {vsitestr},
1361 "Convert atoms to virtual sites" },
1362 { "-heavyh", FALSE0, etBOOL, {&bHeavyH},
1363 "Make hydrogen atoms heavy" },
1364 { "-deuterate", FALSE0, etBOOL, {&bDeuterate},
1365 "Change the mass of hydrogens to 2 amu" },
1366 { "-chargegrp", TRUE1, etBOOL, {&bChargeGroups},
1367 "Use charge groups in the [TT].rtp[tt] file" },
1368 { "-cmap", TRUE1, etBOOL, {&bCmap},
1369 "Use cmap torsions (if enabled in the [TT].rtp[tt] file)" },
1370 { "-renum", TRUE1, etBOOL, {&bRenumRes},
1371 "Renumber the residues consecutively in the output" },
1372 { "-rtpres", TRUE1, etBOOL, {&bRTPresname},
1373 "Use [TT].rtp[tt] entry names as residue names" }
1374 };
1375#define NPARGS((int)(sizeof(pa)/sizeof((pa)[0]))) asize(pa)((int)(sizeof(pa)/sizeof((pa)[0])))
1376
1377 if (!parse_common_args(&argc, argv, 0, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc,
1
Taking false branch
→
1378 0, NULL((void*)0), &oenv))
1379 {
1380 return 0;
1381 }
1382
1383 /* Force field selection, interactive or direct */
1384 choose_ff(strcmp(ff, "select")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(ff) && __builtin_constant_p ("select") && (
__s1_len = strlen (ff), __s2_len = strlen ("select"), (!((size_t
)(const void *)((ff) + 1) - (size_t)(const void *)(ff) == 1) ||
__s1_len >= 4) && (!((size_t)(const void *)(("select"
) + 1) - (size_t)(const void *)("select") == 1) || __s2_len >=
4)) ? __builtin_strcmp (ff, "select") : (__builtin_constant_p
(ff) && ((size_t)(const void *)((ff) + 1) - (size_t)
(const void *)(ff) == 1) && (__s1_len = strlen (ff), __s1_len
< 4) ? (__builtin_constant_p ("select") && ((size_t
)(const void *)(("select") + 1) - (size_t)(const void *)("select"
) == 1) ? __builtin_strcmp (ff, "select") : (__extension__ ({
const unsigned char *__s2 = (const unsigned char *) (const char
*) ("select"); int __result = (((const unsigned char *) (const
char *) (ff))[0] - __s2[0]); if (__s1_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
ff))[1] - __s2[1]); if (__s1_len > 1 && __result ==
0) { __result = (((const unsigned char *) (const char *) (ff
))[2] - __s2[2]); if (__s1_len > 2 && __result == 0
) __result = (((const unsigned char *) (const char *) (ff))[3
] - __s2[3]); } } __result; }))) : (__builtin_constant_p ("select"
) && ((size_t)(const void *)(("select") + 1) - (size_t
)(const void *)("select") == 1) && (__s2_len = strlen
("select"), __s2_len < 4) ? (__builtin_constant_p (ff) &&
((size_t)(const void *)((ff) + 1) - (size_t)(const void *)(ff
) == 1) ? __builtin_strcmp (ff, "select") : (- (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (ff); int __result = (((const unsigned char *) (const
char *) ("select"))[0] - __s2[0]); if (__s2_len > 0 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) ("select"))[1] - __s2[1]); if (__s2_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) ("select"))[2] - __s2[2]); if (__s2_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) ("select"))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(ff, "select")))); }) == 0 ? NULL((void*)0) : ff,
1385 forcefield, sizeof(forcefield),
1386 ffdir, sizeof(ffdir));
1387
1388 if (strlen(forcefield) > 0)
2
←
Taking true branch
→
1389 {
1390 strcpy(ffname, forcefield);
1391 ffname[0] = toupper(ffname[0])(__extension__ ({ int __res; if (sizeof (ffname[0]) > 1) {
if (__builtin_constant_p (ffname[0])) { int __c = (ffname[0]
); __res = __c < -128 || __c > 255 ? __c : (*__ctype_toupper_loc
())[__c]; } else __res = toupper (ffname[0]); } else __res =
(*__ctype_toupper_loc ())[(int) (ffname[0])]; __res; }));
1392 }
1393 else
1394 {
1395 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1395, "Empty forcefield string");
1396 }
1397
1398 printf("\nUsing the %s force field in directory %s\n\n",
1399 ffname, ffdir);
1400
1401 choose_watermodel(watstr[0], ffdir, &watermodel);
1402
1403 if (bInter)
3
←
Taking false branch
→
1404 {
1405 /* if anything changes here, also change description of -inter */
1406 bCysMan = TRUE1;
1407 bTerMan = TRUE1;
1408 bLysMan = TRUE1;
1409 bArgMan = TRUE1;
1410 bAspMan = TRUE1;
1411 bGluMan = TRUE1;
1412 bGlnMan = TRUE1;
1413 bHisMan = TRUE1;
1414 }
1415
1416 if (bHeavyH)
4
←
Taking false branch
→
1417 {
1418 mHmult = 4.0;
1419 }
1420 else if (bDeuterate)
5
←
Taking false branch
→
1421 {
1422 mHmult = 2.0;
1423 }
1424 else
1425 {
1426 mHmult = 1.0;
1427 }
1428
1429 switch (vsitestr[0][0])
6
←
Array access results in a null pointer dereference
1430 {
1431 case 'n': /* none */
1432 bVsites = FALSE0;
1433 bVsiteAromatics = FALSE0;
1434 break;
1435 case 'h': /* hydrogens */
1436 bVsites = TRUE1;
1437 bVsiteAromatics = FALSE0;
1438 break;
1439 case 'a': /* aromatics */
1440 bVsites = TRUE1;
1441 bVsiteAromatics = TRUE1;
1442 break;
1443 default:
1444 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1444, "DEATH HORROR in $s (%d): vsitestr[0]='%s'",
1445 __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c", __LINE__1445, vsitestr[0]);
1446 } /* end switch */
1447
1448 /* Open the symbol table */
1449 open_symtab(&symtab);
1450
1451 /* Residue type database */
1452 gmx_residuetype_init(&rt);
1453
1454 /* Read residue renaming database(s), if present */
1455 nrrn = fflib_search_file_end(ffdir, ".r2b", FALSE0, &rrn);
1456
1457 nrtprename = 0;
1458 rtprename = NULL((void*)0);
1459 for (i = 0; i < nrrn; i++)
1460 {
1461 fp = fflib_open(rrn[i]);
1462 read_rtprename(rrn[i], fp, &nrtprename, &rtprename);
1463 gmx_ffclose(fp);
1464 sfree(rrn[i])save_free("rrn[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1464, (rrn[i]));
1465 }
1466 sfree(rrn)save_free("rrn", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1466, (rrn));
1467
1468 /* Add all alternative names from the residue renaming database to the list of recognized amino/nucleic acids. */
1469 naa = 0;
1470 for (i = 0; i < nrtprename; i++)
1471 {
1472 rc = gmx_residuetype_get_type(rt, rtprename[i].gmx, &p_restype);
1473
1474 /* Only add names if the 'standard' gromacs/iupac base name was found */
1475 if (rc == 0)
1476 {
1477 gmx_residuetype_add(rt, rtprename[i].main, p_restype);
1478 gmx_residuetype_add(rt, rtprename[i].nter, p_restype);
1479 gmx_residuetype_add(rt, rtprename[i].cter, p_restype);
1480 gmx_residuetype_add(rt, rtprename[i].bter, p_restype);
1481 }
1482 }
1483
1484 clear_mat(box);
1485 if (watermodel != NULL((void*)0) && (strstr(watermodel, "4p") ||
1486 strstr(watermodel, "4P")))
1487 {
1488 watres = "HO4";
1489 }
1490 else if (watermodel != NULL((void*)0) && (strstr(watermodel, "5p") ||
1491 strstr(watermodel, "5P")))
1492 {
1493 watres = "HO5";
1494 }
1495 else
1496 {
1497 watres = "HOH";
1498 }
1499
1500 aps = gmx_atomprop_init();
1501 natom = read_pdball(opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn_null("-q", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), title,
1502 &pdba_all, &pdbx, &ePBC, box, bRemoveH, &symtab, rt, watres,
1503 aps, bVerbose);
1504
1505 if (natom == 0)
1506 {
1507 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1507, "No atoms found in pdb file %s\n", opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm));
1508 }
1509
1510 printf("Analyzing pdb file\n");
1511 nch = 0;
1512 maxch = 0;
1513 nwaterchain = 0;
1514
1515 modify_chain_numbers(&pdba_all, chainsep[0]);
1516
1517 nchainmerges = 0;
1518
1519 this_atomname = NULL((void*)0);
1520 this_atomnum = -1;
1521 this_resname = NULL((void*)0);
1522 this_resnum = -1;
1523 this_chainid = '?';
1524 this_chainnumber = -1;
1525 this_chainstart = 0;
1526 /* Keep the compiler happy */
1527 prev_chainstart = 0;
1528
1529 pdb_ch = NULL((void*)0);
1530
1531 bMerged = FALSE0;
1532 for (i = 0; (i < natom); i++)
1533 {
1534 ri = &pdba_all.resinfo[pdba_all.atom[i].resind];
1535
1536 prev_atomname = this_atomname;
1537 prev_atomnum = this_atomnum;
1538 prev_resname = this_resname;
1539 prev_resnum = this_resnum;
1540 prev_chainid = this_chainid;
1541 prev_chainnumber = this_chainnumber;
1542 if (!bMerged)
1543 {
1544 prev_chainstart = this_chainstart;
1545 }
1546
1547 this_atomname = *pdba_all.atomname[i];
1548 this_atomnum = (pdba_all.pdbinfo != NULL((void*)0)) ? pdba_all.pdbinfo[i].atomnr : i+1;
1549 this_resname = *ri->name;
1550 this_resnum = ri->nr;
1551 this_chainid = ri->chainid;
1552 this_chainnumber = ri->chainnum;
1553
1554 bWat = gmx_strcasecmp(*ri->name, watres) == 0;
1555 if ((i == 0) || (this_chainnumber != prev_chainnumber) || (bWat != bPrevWat))
1556 {
1557 this_chainstart = pdba_all.atom[i].resind;
1558
1559 bMerged = FALSE0;
1560 if (i > 0 && !bWat)
1561 {
1562 if (!strncmp(merge[0], "int", 3)(__extension__ (__builtin_constant_p (3) && ((__builtin_constant_p
(merge[0]) && strlen (merge[0]) < ((size_t) (3)))
|| (__builtin_constant_p ("int") && strlen ("int") <
((size_t) (3)))) ? __extension__ ({ size_t __s1_len, __s2_len
; (__builtin_constant_p (merge[0]) && __builtin_constant_p
("int") && (__s1_len = strlen (merge[0]), __s2_len =
strlen ("int"), (!((size_t)(const void *)((merge[0]) + 1) - (
size_t)(const void *)(merge[0]) == 1) || __s1_len >= 4) &&
(!((size_t)(const void *)(("int") + 1) - (size_t)(const void
*)("int") == 1) || __s2_len >= 4)) ? __builtin_strcmp (merge
[0], "int") : (__builtin_constant_p (merge[0]) && ((size_t
)(const void *)((merge[0]) + 1) - (size_t)(const void *)(merge
[0]) == 1) && (__s1_len = strlen (merge[0]), __s1_len
< 4) ? (__builtin_constant_p ("int") && ((size_t)
(const void *)(("int") + 1) - (size_t)(const void *)("int") ==
1) ? __builtin_strcmp (merge[0], "int") : (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
("int"); int __result = (((const unsigned char *) (const char
*) (merge[0]))[0] - __s2[0]); if (__s1_len > 0 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (merge[0]))[1] - __s2[1]); if (__s1_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (merge[0]))[2] - __s2[2]); if (__s1_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) (merge[0]))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
("int") && ((size_t)(const void *)(("int") + 1) - (size_t
)(const void *)("int") == 1) && (__s2_len = strlen ("int"
), __s2_len < 4) ? (__builtin_constant_p (merge[0]) &&
((size_t)(const void *)((merge[0]) + 1) - (size_t)(const void
*)(merge[0]) == 1) ? __builtin_strcmp (merge[0], "int") : (-
(__extension__ ({ const unsigned char *__s2 = (const unsigned
char *) (const char *) (merge[0]); int __result = (((const unsigned
char *) (const char *) ("int"))[0] - __s2[0]); if (__s2_len >
0 && __result == 0) { __result = (((const unsigned char
*) (const char *) ("int"))[1] - __s2[1]); if (__s2_len > 1
&& __result == 0) { __result = (((const unsigned char
*) (const char *) ("int"))[2] - __s2[2]); if (__s2_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) ("int"))[3] - __s2[3]); } } __result; })))) :
__builtin_strcmp (merge[0], "int")))); }) : strncmp (merge[0
], "int", 3))))
1563 {
1564 printf("Merge chain ending with residue %s%d (chain id '%c', atom %d %s) and chain starting with\n"
1565 "residue %s%d (chain id '%c', atom %d %s) into a single moleculetype (keeping termini)? [n/y]\n",
1566 prev_resname, prev_resnum, prev_chainid, prev_atomnum, prev_atomname,
1567 this_resname, this_resnum, this_chainid, this_atomnum, this_atomname);
1568
1569 if (NULL((void*)0) == fgets(select, STRLEN4096-1, stdinstdin))
1570 {
1571 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1571, "Error reading from stdin");
1572 }
1573 bMerged = (select[0] == 'y');
1574 }
1575 else if (!strncmp(merge[0], "all", 3)(__extension__ (__builtin_constant_p (3) && ((__builtin_constant_p
(merge[0]) && strlen (merge[0]) < ((size_t) (3)))
|| (__builtin_constant_p ("all") && strlen ("all") <
((size_t) (3)))) ? __extension__ ({ size_t __s1_len, __s2_len
; (__builtin_constant_p (merge[0]) && __builtin_constant_p
("all") && (__s1_len = strlen (merge[0]), __s2_len =
strlen ("all"), (!((size_t)(const void *)((merge[0]) + 1) - (
size_t)(const void *)(merge[0]) == 1) || __s1_len >= 4) &&
(!((size_t)(const void *)(("all") + 1) - (size_t)(const void
*)("all") == 1) || __s2_len >= 4)) ? __builtin_strcmp (merge
[0], "all") : (__builtin_constant_p (merge[0]) && ((size_t
)(const void *)((merge[0]) + 1) - (size_t)(const void *)(merge
[0]) == 1) && (__s1_len = strlen (merge[0]), __s1_len
< 4) ? (__builtin_constant_p ("all") && ((size_t)
(const void *)(("all") + 1) - (size_t)(const void *)("all") ==
1) ? __builtin_strcmp (merge[0], "all") : (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
("all"); int __result = (((const unsigned char *) (const char
*) (merge[0]))[0] - __s2[0]); if (__s1_len > 0 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (merge[0]))[1] - __s2[1]); if (__s1_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (merge[0]))[2] - __s2[2]); if (__s1_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) (merge[0]))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
("all") && ((size_t)(const void *)(("all") + 1) - (size_t
)(const void *)("all") == 1) && (__s2_len = strlen ("all"
), __s2_len < 4) ? (__builtin_constant_p (merge[0]) &&
((size_t)(const void *)((merge[0]) + 1) - (size_t)(const void
*)(merge[0]) == 1) ? __builtin_strcmp (merge[0], "all") : (-
(__extension__ ({ const unsigned char *__s2 = (const unsigned
char *) (const char *) (merge[0]); int __result = (((const unsigned
char *) (const char *) ("all"))[0] - __s2[0]); if (__s2_len >
0 && __result == 0) { __result = (((const unsigned char
*) (const char *) ("all"))[1] - __s2[1]); if (__s2_len > 1
&& __result == 0) { __result = (((const unsigned char
*) (const char *) ("all"))[2] - __s2[2]); if (__s2_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) ("all"))[3] - __s2[3]); } } __result; })))) :
__builtin_strcmp (merge[0], "all")))); }) : strncmp (merge[0
], "all", 3))))
1576 {
1577 bMerged = TRUE1;
1578 }
1579 }
1580
1581 if (bMerged)
1582 {
1583 pdb_ch[nch-1].chainstart[pdb_ch[nch-1].nterpairs] =
1584 pdba_all.atom[i].resind - prev_chainstart;
1585 pdb_ch[nch-1].nterpairs++;
1586 srenew(pdb_ch[nch-1].chainstart, pdb_ch[nch-1].nterpairs+1)(pdb_ch[nch-1].chainstart) = save_realloc("pdb_ch[nch-1].chainstart"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1586, (pdb_ch[nch-1].chainstart), (pdb_ch[nch-1].nterpairs+
1), sizeof(*(pdb_ch[nch-1].chainstart)));
1587 nchainmerges++;
1588 }
1589 else
1590 {
1591 /* set natom for previous chain */
1592 if (nch > 0)
1593 {
1594 pdb_ch[nch-1].natom = i-pdb_ch[nch-1].start;
1595 }
1596 if (bWat)
1597 {
1598 nwaterchain++;
1599 ri->chainid = ' ';
1600 }
1601 /* check if chain identifier was used before */
1602 for (j = 0; (j < nch); j++)
1603 {
1604 if (pdb_ch[j].chainid != ' ' && pdb_ch[j].chainid == ri->chainid)
1605 {
1606 printf("WARNING: Chain identifier '%c' is used in two non-sequential blocks.\n"
1607 "They will be treated as separate chains unless you reorder your file.\n",
1608 ri->chainid);
1609 }
1610 }
1611 if (nch == maxch)
1612 {
1613 maxch += 16;
1614 srenew(pdb_ch, maxch)(pdb_ch) = save_realloc("pdb_ch", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1614, (pdb_ch), (maxch), sizeof(*(pdb_ch)));
1615 }
1616 pdb_ch[nch].chainid = ri->chainid;
1617 pdb_ch[nch].chainnum = ri->chainnum;
1618 pdb_ch[nch].start = i;
1619 pdb_ch[nch].bAllWat = bWat;
1620 if (bWat)
1621 {
1622 pdb_ch[nch].nterpairs = 0;
1623 }
1624 else
1625 {
1626 pdb_ch[nch].nterpairs = 1;
1627 }
1628 snew(pdb_ch[nch].chainstart, pdb_ch[nch].nterpairs+1)(pdb_ch[nch].chainstart) = save_calloc("pdb_ch[nch].chainstart"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1628, (pdb_ch[nch].nterpairs+1), sizeof(*(pdb_ch[nch].chainstart
)));
1629 /* modified [nch] to [0] below */
1630 pdb_ch[nch].chainstart[0] = 0;
1631 nch++;
1632 }
1633 }
1634 bPrevWat = bWat;
1635 }
1636 pdb_ch[nch-1].natom = natom-pdb_ch[nch-1].start;
1637
1638 /* set all the water blocks at the end of the chain */
1639 snew(swap_index, nch)(swap_index) = save_calloc("swap_index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1639, (nch), sizeof(*(swap_index)));
1640 j = 0;
1641 for (i = 0; i < nch; i++)
1642 {
1643 if (!pdb_ch[i].bAllWat)
1644 {
1645 swap_index[j] = i;
1646 j++;
1647 }
1648 }
1649 for (i = 0; i < nch; i++)
1650 {
1651 if (pdb_ch[i].bAllWat)
1652 {
1653 swap_index[j] = i;
1654 j++;
1655 }
1656 }
1657 if (nwaterchain > 1)
1658 {
1659 printf("Moved all the water blocks to the end\n");
1660 }
1661
1662 snew(chains, nch)(chains) = save_calloc("chains", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1662, (nch), sizeof(*(chains)));
1663 /* copy pdb data and x for all chains */
1664 for (i = 0; (i < nch); i++)
1665 {
1666 si = swap_index[i];
1667 chains[i].chainid = pdb_ch[si].chainid;
1668 chains[i].chainnum = pdb_ch[si].chainnum;
1669 chains[i].bAllWat = pdb_ch[si].bAllWat;
1670 chains[i].nterpairs = pdb_ch[si].nterpairs;
1671 chains[i].chainstart = pdb_ch[si].chainstart;
1672 snew(chains[i].ntdb, pdb_ch[si].nterpairs)(chains[i].ntdb) = save_calloc("chains[i].ntdb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1672, (pdb_ch[si].nterpairs), sizeof(*(chains[i].ntdb)));
1673 snew(chains[i].ctdb, pdb_ch[si].nterpairs)(chains[i].ctdb) = save_calloc("chains[i].ctdb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1673, (pdb_ch[si].nterpairs), sizeof(*(chains[i].ctdb)));
1674 snew(chains[i].r_start, pdb_ch[si].nterpairs)(chains[i].r_start) = save_calloc("chains[i].r_start", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1674, (pdb_ch[si].nterpairs), sizeof(*(chains[i].r_start)));
1675 snew(chains[i].r_end, pdb_ch[si].nterpairs)(chains[i].r_end) = save_calloc("chains[i].r_end", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1675, (pdb_ch[si].nterpairs), sizeof(*(chains[i].r_end)));
1676
1677 snew(chains[i].pdba, 1)(chains[i].pdba) = save_calloc("chains[i].pdba", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1677, (1), sizeof(*(chains[i].pdba)));
1678 init_t_atoms(chains[i].pdba, pdb_ch[si].natom, TRUE1);
1679 snew(chains[i].x, chains[i].pdba->nr)(chains[i].x) = save_calloc("chains[i].x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1679, (chains[i].pdba->nr), sizeof(*(chains[i].x)));
1680 for (j = 0; j < chains[i].pdba->nr; j++)
1681 {
1682 chains[i].pdba->atom[j] = pdba_all.atom[pdb_ch[si].start+j];
1683 snew(chains[i].pdba->atomname[j], 1)(chains[i].pdba->atomname[j]) = save_calloc("chains[i].pdba->atomname[j]"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1683, (1), sizeof(*(chains[i].pdba->atomname[j])));
1684 *chains[i].pdba->atomname[j] =
1685 strdup(*pdba_all.atomname[pdb_ch[si].start+j])(__extension__ (__builtin_constant_p (*pdba_all.atomname[pdb_ch
[si].start+j]) && ((size_t)(const void *)((*pdba_all.
atomname[pdb_ch[si].start+j]) + 1) - (size_t)(const void *)(*
pdba_all.atomname[pdb_ch[si].start+j]) == 1) ? (((const char *
) (*pdba_all.atomname[pdb_ch[si].start+j]))[0] == '\0' ? (char
*) calloc ((size_t) 1, (size_t) 1) : ({ size_t __len = strlen
(*pdba_all.atomname[pdb_ch[si].start+j]) + 1; char *__retval
= (char *) malloc (__len); if (__retval != ((void*)0)) __retval
= (char *) memcpy (__retval, *pdba_all.atomname[pdb_ch[si].start
+j], __len); __retval; })) : __strdup (*pdba_all.atomname[pdb_ch
[si].start+j])));
1686 chains[i].pdba->pdbinfo[j] = pdba_all.pdbinfo[pdb_ch[si].start+j];
1687 copy_rvec(pdbx[pdb_ch[si].start+j], chains[i].x[j]);
1688 }
1689 /* Re-index the residues assuming that the indices are continuous */
1690 k = chains[i].pdba->atom[0].resind;
1691 nres = chains[i].pdba->atom[chains[i].pdba->nr-1].resind - k + 1;
1692 chains[i].pdba->nres = nres;
1693 for (j = 0; j < chains[i].pdba->nr; j++)
1694 {
1695 chains[i].pdba->atom[j].resind -= k;
1696 }
1697 srenew(chains[i].pdba->resinfo, nres)(chains[i].pdba->resinfo) = save_realloc("chains[i].pdba->resinfo"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1697, (chains[i].pdba->resinfo), (nres), sizeof(*(chains
[i].pdba->resinfo)));
1698 for (j = 0; j < nres; j++)
1699 {
1700 chains[i].pdba->resinfo[j] = pdba_all.resinfo[k+j];
1701 snew(chains[i].pdba->resinfo[j].name, 1)(chains[i].pdba->resinfo[j].name) = save_calloc("chains[i].pdba->resinfo[j].name"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1701, (1), sizeof(*(chains[i].pdba->resinfo[j].name)));
1702 *chains[i].pdba->resinfo[j].name = strdup(*pdba_all.resinfo[k+j].name)(__extension__ (__builtin_constant_p (*pdba_all.resinfo[k+j].
name) && ((size_t)(const void *)((*pdba_all.resinfo[k
+j].name) + 1) - (size_t)(const void *)(*pdba_all.resinfo[k+j
].name) == 1) ? (((const char *) (*pdba_all.resinfo[k+j].name
))[0] == '\0' ? (char *) calloc ((size_t) 1, (size_t) 1) : ({
size_t __len = strlen (*pdba_all.resinfo[k+j].name) + 1; char
*__retval = (char *) malloc (__len); if (__retval != ((void*
)0)) __retval = (char *) memcpy (__retval, *pdba_all.resinfo[
k+j].name, __len); __retval; })) : __strdup (*pdba_all.resinfo
[k+j].name)));
1703 /* make all chain identifiers equal to that of the chain */
1704 chains[i].pdba->resinfo[j].chainid = pdb_ch[si].chainid;
1705 }
1706 }
1707
1708 if (nchainmerges > 0)
1709 {
1710 printf("\nMerged chains into joint molecule definitions at %d places.\n\n",
1711 nchainmerges);
1712 }
1713
1714 printf("There are %d chains and %d blocks of water and "
1715 "%d residues with %d atoms\n",
1716 nch-nwaterchain, nwaterchain,
1717 pdba_all.resinfo[pdba_all.atom[natom-1].resind].nr, natom);
1718
1719 printf("\n %5s %4s %6s\n", "chain", "#res", "#atoms");
1720 for (i = 0; (i < nch); i++)
1721 {
1722 printf(" %d '%c' %5d %6d %s\n",
1723 i+1, chains[i].chainid ? chains[i].chainid : '-',
1724 chains[i].pdba->nres, chains[i].pdba->nr,
1725 chains[i].bAllWat ? "(only water)" : "");
1726 }
1727 printf("\n");
1728
1729 check_occupancy(&pdba_all, opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), bVerbose);
1730
1731 /* Read atomtypes... */
1732 atype = read_atype(ffdir, &symtab);
1733
1734 /* read residue database */
1735 printf("Reading residue database... (%s)\n", forcefield);
1736 nrtpf = fflib_search_file_end(ffdir, ".rtp", TRUE1, &rtpf);
1737 nrtp = 0;
1738 restp = NULL((void*)0);
1739 for (i = 0; i < nrtpf; i++)
1740 {
1741 read_resall(rtpf[i], &nrtp, &restp, atype, &symtab, FALSE0);
1742 sfree(rtpf[i])save_free("rtpf[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1742, (rtpf[i]));
1743 }
1744 sfree(rtpf)save_free("rtpf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1744, (rtpf));
1745 if (bNewRTP)
1746 {
1747 /* Not correct with multiple rtp input files with different bonded types */
1748 fp = gmx_fio_fopen("new.rtp", "w");
1749 print_resall(fp, nrtp, restp, atype);
1750 gmx_fio_fclose(fp);
1751 }
1752
1753 /* read hydrogen database */
1754 nah = read_h_db(ffdir, &ah);
1755
1756 /* Read Termini database... */
1757 nNtdb = read_ter_db(ffdir, 'n', &ntdb, atype);
1758 nCtdb = read_ter_db(ffdir, 'c', &ctdb, atype);
1759
1760 top_fn = ftp2fn(efTOP, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
1761 top_file = gmx_fio_fopen(top_fn, "w");
1762
1763 print_top_header(top_file, top_fn, FALSE0, ffdir, mHmult);
1764
1765 nincl = 0;
1766 nmol = 0;
1767 incls = NULL((void*)0);
1768 mols = NULL((void*)0);
1769 nres = 0;
1770 for (chain = 0; (chain < nch); chain++)
1771 {
1772 cc = &(chains[chain]);
1773
1774 /* set pdba, natom and nres to the current chain */
1775 pdba = cc->pdba;
1776 x = cc->x;
1777 natom = cc->pdba->nr;
1778 nres = cc->pdba->nres;
1779
1780 if (cc->chainid && ( cc->chainid != ' ' ) )
1781 {
1782 printf("Processing chain %d '%c' (%d atoms, %d residues)\n",
1783 chain+1, cc->chainid, natom, nres);
1784 }
1785 else
1786 {
1787 printf("Processing chain %d (%d atoms, %d residues)\n",
1788 chain+1, natom, nres);
1789 }
1790
1791 process_chain(pdba, x, bUnA, bUnA, bUnA, bLysMan, bAspMan, bGluMan,
1792 bHisMan, bArgMan, bGlnMan, angle, distance, &symtab,
1793 nrtprename, rtprename);
1794
1795 cc->chainstart[cc->nterpairs] = pdba->nres;
1796 j = 0;
1797 for (i = 0; i < cc->nterpairs; i++)
1798 {
1799 find_nc_ter(pdba, cc->chainstart[i], cc->chainstart[i+1],
1800 &(cc->r_start[j]), &(cc->r_end[j]), rt);
1801
1802 if (cc->r_start[j] >= 0 && cc->r_end[j] >= 0)
1803 {
1804 j++;
1805 }
1806 }
1807 cc->nterpairs = j;
1808 if (cc->nterpairs == 0)
1809 {
1810 printf("Problem with chain definition, or missing terminal residues.\n"
1811 "This chain does not appear to contain a recognized chain molecule.\n"
1812 "If this is incorrect, you can edit residuetypes.dat to modify the behavior.\n");
1813 }
1814
1815 /* Check for disulfides and other special bonds */
1816 nssbonds = mk_specbonds(pdba, x, bCysMan, &ssbonds, bVerbose);
1817
1818 if (nrtprename > 0)
1819 {
1820 rename_resrtp(pdba, cc->nterpairs, cc->r_start, cc->r_end, nrtprename, rtprename,
1821 &symtab, bVerbose);
1822 }
1823
1824 if (debug)
1825 {
1826 if (nch == 1)
1827 {
1828 sprintf(fn, "chain.pdb");
1829 }
1830 else
1831 {
1832 sprintf(fn, "chain_%c%d.pdb", cc->chainid, cc->chainnum);
1833 }
1834 write_sto_conf(fn, title, pdba, x, NULL((void*)0), ePBC, box);
1835 }
1836
1837
1838 for (i = 0; i < cc->nterpairs; i++)
1839 {
1840
1841 /* Set termini.
1842 * We first apply a filter so we only have the
1843 * termini that can be applied to the residue in question
1844 * (or a generic terminus if no-residue specific is available).
1845 */
1846 /* First the N terminus */
1847 if (nNtdb > 0)
1848 {
1849 tdblist = filter_ter(nrtp, restp, nNtdb, ntdb,
1850 *pdba->resinfo[cc->r_start[i]].name,
1851 *pdba->resinfo[cc->r_start[i]].rtp,
1852 &ntdblist);
1853 if (ntdblist == 0)
1854 {
1855 printf("No suitable end (N or 5') terminus found in database - assuming this residue\n"
1856 "is already in a terminus-specific form and skipping terminus selection.\n");
1857 cc->ntdb[i] = NULL((void*)0);
1858 }
1859 else
1860 {
1861 if (bTerMan && ntdblist > 1)
1862 {
1863 sprintf(select, "Select start terminus type for %s-%d",
1864 *pdba->resinfo[cc->r_start[i]].name,
1865 pdba->resinfo[cc->r_start[i]].nr);
1866 cc->ntdb[i] = choose_ter(ntdblist, tdblist, select);
1867 }
1868 else
1869 {
1870 cc->ntdb[i] = tdblist[0];
1871 }
1872
1873 printf("Start terminus %s-%d: %s\n",
1874 *pdba->resinfo[cc->r_start[i]].name,
1875 pdba->resinfo[cc->r_start[i]].nr,
1876 (cc->ntdb[i])->name);
1877 sfree(tdblist)save_free("tdblist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1877, (tdblist));
1878 }
1879 }
1880 else
1881 {
1882 cc->ntdb[i] = NULL((void*)0);
1883 }
1884
1885 /* And the C terminus */
1886 if (nCtdb > 0)
1887 {
1888 tdblist = filter_ter(nrtp, restp, nCtdb, ctdb,
1889 *pdba->resinfo[cc->r_end[i]].name,
1890 *pdba->resinfo[cc->r_end[i]].rtp,
1891 &ntdblist);
1892 if (ntdblist == 0)
1893 {
1894 printf("No suitable end (C or 3') terminus found in database - assuming this residue\n"
1895 "is already in a terminus-specific form and skipping terminus selection.\n");
1896 cc->ctdb[i] = NULL((void*)0);
1897 }
1898 else
1899 {
1900 if (bTerMan && ntdblist > 1)
1901 {
1902 sprintf(select, "Select end terminus type for %s-%d",
1903 *pdba->resinfo[cc->r_end[i]].name,
1904 pdba->resinfo[cc->r_end[i]].nr);
1905 cc->ctdb[i] = choose_ter(ntdblist, tdblist, select);
1906 }
1907 else
1908 {
1909 cc->ctdb[i] = tdblist[0];
1910 }
1911 printf("End terminus %s-%d: %s\n",
1912 *pdba->resinfo[cc->r_end[i]].name,
1913 pdba->resinfo[cc->r_end[i]].nr,
1914 (cc->ctdb[i])->name);
1915 sfree(tdblist)save_free("tdblist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1915, (tdblist));
1916 }
1917 }
1918 else
1919 {
1920 cc->ctdb[i] = NULL((void*)0);
1921 }
1922 }
1923 /* lookup hackblocks and rtp for all residues */
1924 get_hackblocks_rtp(&hb_chain, &restp_chain,
1925 nrtp, restp, pdba->nres, pdba->resinfo,
1926 cc->nterpairs, cc->ntdb, cc->ctdb, cc->r_start, cc->r_end);
1927 /* ideally, now we would not need the rtp itself anymore, but do
1928 everything using the hb and restp arrays. Unfortunately, that
1929 requires some re-thinking of code in gen_vsite.c, which I won't
1930 do now :( AF 26-7-99 */
1931
1932 rename_atoms(NULL((void*)0), ffdir,
1933 pdba, &symtab, restp_chain, FALSE0, rt, FALSE0, bVerbose);
1934
1935 match_atomnames_with_rtp(restp_chain, hb_chain, pdba, x, bVerbose);
1936
1937 if (bSort)
1938 {
1939 block = new_blocka();
1940 snew(gnames, 1)(gnames) = save_calloc("gnames", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1940, (1), sizeof(*(gnames)));
1941 sort_pdbatoms(restp_chain, natom, &pdba, &x, block, &gnames);
1942 natom = remove_duplicate_atoms(pdba, x, bVerbose);
1943 if (ftp2bSet(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
1944 {
1945 if (bRemoveH)
1946 {
1947 fprintf(stderrstderr, "WARNING: with the -remh option the generated "
1948 "index file (%s) might be useless\n"
1949 "(the index file is generated before hydrogens are added)",
1950 ftp2fn(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm));
1951 }
1952 write_index(ftp2fn(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), block, gnames, FALSE0, 0);
1953 }
1954 for (i = 0; i < block->nr; i++)
1955 {
1956 sfree(gnames[i])save_free("gnames[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1956, (gnames[i]));
1957 }
1958 sfree(gnames)save_free("gnames", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 1958, (gnames));
1959 done_blocka(block);
1960 }
1961 else
1962 {
1963 fprintf(stderrstderr, "WARNING: "
1964 "without sorting no check for duplicate atoms can be done\n");
1965 }
1966
1967 /* Generate Hydrogen atoms (and termini) in the sequence */
1968 printf("Generating any missing hydrogen atoms and/or adding termini.\n");
1969 natom = add_h(&pdba, &x, nah, ah,
1970 cc->nterpairs, cc->ntdb, cc->ctdb, cc->r_start, cc->r_end, bAllowMissing,
1971 NULL((void*)0), NULL((void*)0), TRUE1, FALSE0);
1972 printf("Now there are %d residues with %d atoms\n",
1973 pdba->nres, pdba->nr);
1974 if (debug)
1975 {
1976 write_pdbfile(debug, title, pdba, x, ePBC, box, ' ', 0, NULL((void*)0), TRUE1);
1977 }
1978
1979 if (debug)
1980 {
1981 for (i = 0; (i < natom); i++)
1982 {
1983 fprintf(debug, "Res %s%d atom %d %s\n",
1984 *(pdba->resinfo[pdba->atom[i].resind].name),
1985 pdba->resinfo[pdba->atom[i].resind].nr, i+1, *pdba->atomname[i]);
1986 }
1987 }
1988
1989 strcpy(posre_fn, ftp2fn(efITP, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm));
1990
1991 /* make up molecule name(s) */
1992
1993 k = (cc->nterpairs > 0 && cc->r_start[0] >= 0) ? cc->r_start[0] : 0;
1994
1995 gmx_residuetype_get_type(rt, *pdba->resinfo[k].name, &p_restype);
1996
1997 suffix[0] = '\0';
1998
1999 if (cc->bAllWat)
2000 {
2001 sprintf(molname, "Water");
2002 }
2003 else
2004 {
2005 this_chainid = cc->chainid;
2006
2007 /* Add the chain id if we have one */
2008 if (this_chainid != ' ')
2009 {
2010 sprintf(buf, "_chain_%c", this_chainid);
2011 strcat(suffix, buf);
2012 }
2013
2014 /* Check if there have been previous chains with the same id */
2015 nid_used = 0;
2016 for (k = 0; k < chain; k++)
2017 {
2018 if (cc->chainid == chains[k].chainid)
2019 {
2020 nid_used++;
2021 }
2022 }
2023 /* Add the number for this chain identifier if there are multiple copies */
2024 if (nid_used > 0)
2025 {
2026
2027 sprintf(buf, "%d", nid_used+1);
2028 strcat(suffix, buf);
2029 }
2030
2031 if (strlen(suffix) > 0)
2032 {
2033 sprintf(molname, "%s%s", p_restype, suffix);
2034 }
2035 else
2036 {
2037 strcpy(molname, p_restype);
2038 }
2039 }
2040
2041 if ((nch-nwaterchain > 1) && !cc->bAllWat)
2042 {
2043 bITP = TRUE1;
2044 strcpy(itp_fn, top_fn);
2045 printf("Chain time...\n");
2046 c = strrchr(itp_fn, '.');
2047 sprintf(c, "_%s.itp", molname);
2048 c = strrchr(posre_fn, '.');
2049 sprintf(c, "_%s.itp", molname);
2050 if (strcmp(itp_fn, posre_fn)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(itp_fn) && __builtin_constant_p (posre_fn) &&
(__s1_len = strlen (itp_fn), __s2_len = strlen (posre_fn), (
!((size_t)(const void *)((itp_fn) + 1) - (size_t)(const void *
)(itp_fn) == 1) || __s1_len >= 4) && (!((size_t)(const
void *)((posre_fn) + 1) - (size_t)(const void *)(posre_fn) ==
1) || __s2_len >= 4)) ? __builtin_strcmp (itp_fn, posre_fn
) : (__builtin_constant_p (itp_fn) && ((size_t)(const
void *)((itp_fn) + 1) - (size_t)(const void *)(itp_fn) == 1)
&& (__s1_len = strlen (itp_fn), __s1_len < 4) ? (
__builtin_constant_p (posre_fn) && ((size_t)(const void
*)((posre_fn) + 1) - (size_t)(const void *)(posre_fn) == 1) ?
__builtin_strcmp (itp_fn, posre_fn) : (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
(posre_fn); int __result = (((const unsigned char *) (const char
*) (itp_fn))[0] - __s2[0]); if (__s1_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
itp_fn))[1] - __s2[1]); if (__s1_len > 1 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
itp_fn))[2] - __s2[2]); if (__s1_len > 2 && __result
== 0) __result = (((const unsigned char *) (const char *) (itp_fn
))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p (
posre_fn) && ((size_t)(const void *)((posre_fn) + 1) -
(size_t)(const void *)(posre_fn) == 1) && (__s2_len =
strlen (posre_fn), __s2_len < 4) ? (__builtin_constant_p (
itp_fn) && ((size_t)(const void *)((itp_fn) + 1) - (size_t
)(const void *)(itp_fn) == 1) ? __builtin_strcmp (itp_fn, posre_fn
) : (- (__extension__ ({ const unsigned char *__s2 = (const unsigned
char *) (const char *) (itp_fn); int __result = (((const unsigned
char *) (const char *) (posre_fn))[0] - __s2[0]); if (__s2_len
> 0 && __result == 0) { __result = (((const unsigned
char *) (const char *) (posre_fn))[1] - __s2[1]); if (__s2_len
> 1 && __result == 0) { __result = (((const unsigned
char *) (const char *) (posre_fn))[2] - __s2[2]); if (__s2_len
> 2 && __result == 0) __result = (((const unsigned
char *) (const char *) (posre_fn))[3] - __s2[3]); } } __result
; })))) : __builtin_strcmp (itp_fn, posre_fn)))); }) == 0)
2051 {
2052 strcpy(buf_fn, posre_fn);
2053 c = strrchr(buf_fn, '.');
2054 *c = '\0';
2055 sprintf(posre_fn, "%s_pr.itp", buf_fn);
2056 }
2057
2058 nincl++;
2059 srenew(incls, nincl)(incls) = save_realloc("incls", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 2059, (incls), (nincl), sizeof(*(incls)));
2060 incls[nincl-1] = strdup(itp_fn)(__extension__ (__builtin_constant_p (itp_fn) && ((size_t
)(const void *)((itp_fn) + 1) - (size_t)(const void *)(itp_fn
) == 1) ? (((const char *) (itp_fn))[0] == '\0' ? (char *) calloc
((size_t) 1, (size_t) 1) : ({ size_t __len = strlen (itp_fn)
+ 1; char *__retval = (char *) malloc (__len); if (__retval !=
((void*)0)) __retval = (char *) memcpy (__retval, itp_fn, __len
); __retval; })) : __strdup (itp_fn)));
2061 itp_file = gmx_fio_fopen(itp_fn, "w");
2062 }
2063 else
2064 {
2065 bITP = FALSE0;
2066 }
2067
2068 srenew(mols, nmol+1)(mols) = save_realloc("mols", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 2068, (mols), (nmol+1), sizeof(*(mols)));
2069 if (cc->bAllWat)
2070 {
2071 mols[nmol].name = strdup("SOL")(__extension__ (__builtin_constant_p ("SOL") && ((size_t
)(const void *)(("SOL") + 1) - (size_t)(const void *)("SOL") ==
1) ? (((const char *) ("SOL"))[0] == '\0' ? (char *) calloc (
(size_t) 1, (size_t) 1) : ({ size_t __len = strlen ("SOL") + 1
; char *__retval = (char *) malloc (__len); if (__retval != (
(void*)0)) __retval = (char *) memcpy (__retval, "SOL", __len
); __retval; })) : __strdup ("SOL")));
2072 mols[nmol].nr = pdba->nres;
2073 }
2074 else
2075 {
2076 mols[nmol].name = strdup(molname)(__extension__ (__builtin_constant_p (molname) && ((size_t
)(const void *)((molname) + 1) - (size_t)(const void *)(molname
) == 1) ? (((const char *) (molname))[0] == '\0' ? (char *) calloc
((size_t) 1, (size_t) 1) : ({ size_t __len = strlen (molname
) + 1; char *__retval = (char *) malloc (__len); if (__retval
!= ((void*)0)) __retval = (char *) memcpy (__retval, molname
, __len); __retval; })) : __strdup (molname)));
2077 mols[nmol].nr = 1;
2078 }
2079 nmol++;
2080
2081 if (bITP)
2082 {
2083 print_top_comment(itp_file, itp_fn, ffdir, TRUE1);
2084 }
2085
2086 if (cc->bAllWat)
2087 {
2088 top_file2 = NULL((void*)0);
2089 }
2090 else
2091 if (bITP)
2092 {
2093 top_file2 = itp_file;
2094 }
2095 else
2096 {
2097 top_file2 = top_file;
2098 }
2099
2100 pdb2top(top_file2, posre_fn, molname, pdba, &x, atype, &symtab,
2101 nrtp, restp,
2102 restp_chain, hb_chain,
2103 bAllowMissing,
2104 bVsites, bVsiteAromatics, ffdir,
2105 mHmult, nssbonds, ssbonds,
2106 long_bond_dist, short_bond_dist, bDeuterate, bChargeGroups, bCmap,
2107 bRenumRes, bRTPresname);
2108
2109 if (!cc->bAllWat)
2110 {
2111 write_posres(posre_fn, pdba, posre_fc);
2112 }
2113
2114 if (bITP)
2115 {
2116 gmx_fio_fclose(itp_file);
2117 }
2118
2119 /* pdba and natom have been reassigned somewhere so: */
2120 cc->pdba = pdba;
2121 cc->x = x;
2122
2123 if (debug)
2124 {
2125 if (cc->chainid == ' ')
2126 {
2127 sprintf(fn, "chain.pdb");
2128 }
2129 else
2130 {
2131 sprintf(fn, "chain_%c.pdb", cc->chainid);
2132 }
2133 cool_quote(quote, 255, NULL((void*)0));
2134 write_sto_conf(fn, quote, pdba, x, NULL((void*)0), ePBC, box);
2135 }
2136 }
2137
2138 if (watermodel == NULL((void*)0))
2139 {
2140 for (chain = 0; chain < nch; chain++)
2141 {
2142 if (chains[chain].bAllWat)
2143 {
2144 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 2144, "You have chosen not to include a water model, but there is water in the input file. Select a water model or remove the water from your input file.");
2145 }
2146 }
2147 }
2148 else
2149 {
2150 sprintf(buf_fn, "%s%c%s.itp", ffdir, DIR_SEPARATOR'/', watermodel);
2151 if (!fflib_fexist(buf_fn))
2152 {
2153 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 2153, "The topology file '%s' for the selected water model '%s' can not be found in the force field directory. Select a different water model.",
2154 buf_fn, watermodel);
2155 }
2156 }
2157
2158 print_top_mols(top_file, title, ffdir, watermodel, nincl, incls, nmol, mols);
2159 gmx_fio_fclose(top_file);
2160
2161 gmx_residuetype_destroy(rt);
2162
2163 /* now merge all chains back together */
2164 natom = 0;
2165 nres = 0;
2166 for (i = 0; (i < nch); i++)
2167 {
2168 natom += chains[i].pdba->nr;
2169 nres += chains[i].pdba->nres;
2170 }
2171 snew(atoms, 1)(atoms) = save_calloc("atoms", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 2171, (1), sizeof(*(atoms)));
2172 init_t_atoms(atoms, natom, FALSE0);
2173 for (i = 0; i < atoms->nres; i++)
2174 {
2175 sfree(atoms->resinfo[i].name)save_free("atoms->resinfo[i].name", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 2175, (atoms->resinfo[i].name));
2176 }
2177 sfree(atoms->resinfo)save_free("atoms->resinfo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 2177, (atoms->resinfo));
2178 atoms->nres = nres;
2179 snew(atoms->resinfo, nres)(atoms->resinfo) = save_calloc("atoms->resinfo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 2179, (nres), sizeof(*(atoms->resinfo)));
2180 snew(x, natom)(x) = save_calloc("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/pdb2gmx.c"
, 2180, (natom), sizeof(*(x)));
2181 k = 0;
2182 l = 0;
2183 for (i = 0; (i < nch); i++)
2184 {
2185 if (nch > 1)
2186 {
2187 printf("Including chain %d in system: %d atoms %d residues\n",
2188 i+1, chains[i].pdba->nr, chains[i].pdba->nres);
2189 }
2190 for (j = 0; (j < chains[i].pdba->nr); j++)
2191 {
2192 atoms->atom[k] = chains[i].pdba->atom[j];
2193 atoms->atom[k].resind += l; /* l is processed nr of residues */
2194 atoms->atomname[k] = chains[i].pdba->atomname[j];
2195 atoms->resinfo[atoms->atom[k].resind].chainid = chains[i].chainid;
2196 copy_rvec(chains[i].x[j], x[k]);
2197 k++;
2198 }
2199 for (j = 0; (j < chains[i].pdba->nres); j++)
2200 {
2201 atoms->resinfo[l] = chains[i].pdba->resinfo[j];
2202 if (bRTPresname)
2203 {
2204 atoms->resinfo[l].name = atoms->resinfo[l].rtp;
2205 }
2206 l++;
2207 }
2208 }
2209
2210 if (nch > 1)
2211 {
2212 fprintf(stderrstderr, "Now there are %d atoms and %d residues\n", k, l);
2213 print_sums(atoms, TRUE1);
2214 }
2215
2216 fprintf(stderrstderr, "\nWriting coordinate file...\n");
2217 clear_rvec(box_space);
2218 if (box[0][0] == 0)
2219 {
2220 make_new_box(atoms->nr, x, box, box_space, FALSE0);
2221 }
2222 write_sto_conf(ftp2fn(efSTO, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), title, atoms, x, NULL((void*)0), ePBC, box);
2223
2224 printf("\t\t--------- PLEASE NOTE ------------\n");
2225 printf("You have successfully generated a topology from: %s.\n",
2226 opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm));
2227 if (watermodel != NULL((void*)0))
2228 {
2229 printf("The %s force field and the %s water model are used.\n",
2230 ffname, watermodel);
2231 }
2232 else
2233 {
2234 printf("The %s force field is used.\n",
2235 ffname);
2236 }
2237 printf("\t\t--------- ETON ESAELP ------------\n");
2238
2239 return 0;
2240}