File: | gromacs/gmxpreprocess/pdb2gmx.c |
Location: | line 1769, column 5 |
Description: | Value stored to 'nres' is never read |
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; |
Value stored to 'nres' is never read | |
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 | } |