File: | gromacs/gmxpreprocess/pdb2gmx.c |
Location: | line 823, column 37 |
Description: | Array access (via field 'pdbinfo') results in a null pointer dereference |
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 | ||||
82 | typedef 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 | ||||
91 | static 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 | ||||
108 | static 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 | ||||
131 | static 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 | ||||
145 | static 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 | ||||
159 | static 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 | ||||
173 | static 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 | ||||
187 | static 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 | ||||
201 | static 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 | ||||
213 | static 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 | ||||
258 | static 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 | ||||
307 | static 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 | ||||
362 | static 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 | ||||
375 | static 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 | ||||
393 | static 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 | ||||
412 | static 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 | ||||
435 | static 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 | ||||
490 | void 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 | ||||
515 | static 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 | ||||
580 | void 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 */ | |||
651 | typedef 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 | ||||
659 | int 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 | ||||
684 | static 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 | ||||
776 | static 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 | ||||
839 | void 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 | ||||
909 | static void | |||
910 | modify_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 | ||||
1074 | typedef 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 | ||||
1084 | typedef 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 | ||||
1098 | int 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, | |||
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) | |||
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) | |||
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) | |||
1417 | { | |||
1418 | mHmult = 4.0; | |||
1419 | } | |||
1420 | else if (bDeuterate) | |||
1421 | { | |||
1422 | mHmult = 2.0; | |||
1423 | } | |||
1424 | else | |||
1425 | { | |||
1426 | mHmult = 1.0; | |||
1427 | } | |||
1428 | ||||
1429 | switch (vsitestr[0][0]) | |||
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 | } |