2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017, 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.
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.
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.
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.
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.
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.
48 #include "gromacs/fileio/pdbio.h"
49 #include "gromacs/gmxpreprocess/pdb2top.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/utility/cstringutil.h"
52 #include "gromacs/utility/fatalerror.h"
53 #include "gromacs/utility/smalloc.h"
54 #include "gromacs/utility/strdb.h"
62 c
= toupper(fgetc(stdin
));
64 while ((c
!= 'Y') && (c
!= 'N'));
69 t_specbond
*get_specbonds(int *nspecbond
)
71 const char *sbfile
= "specbond.dat";
73 t_specbond
*sb
= nullptr;
74 char r1buf
[32], r2buf
[32], a1buf
[32], a2buf
[32], nr1buf
[32], nr2buf
[32];
80 nlines
= get_lines(sbfile
, &lines
);
87 for (i
= 0; (i
< nlines
); i
++)
89 if (sscanf(lines
[i
], "%s%s%d%s%s%d%lf%s%s",
90 r1buf
, a1buf
, &nb1
, r2buf
, a2buf
, &nb2
, &length
, nr1buf
, nr2buf
) != 9)
92 fprintf(stderr
, "Invalid line '%s' in %s\n", lines
[i
], sbfile
);
96 sb
[n
].res1
= gmx_strdup(r1buf
);
97 sb
[n
].res2
= gmx_strdup(r2buf
);
98 sb
[n
].newres1
= gmx_strdup(nr1buf
);
99 sb
[n
].newres2
= gmx_strdup(nr2buf
);
100 sb
[n
].atom1
= gmx_strdup(a1buf
);
101 sb
[n
].atom2
= gmx_strdup(a2buf
);
104 sb
[n
].length
= length
;
113 fprintf(stderr
, "%d out of %d lines of %s converted successfully\n",
121 void done_specbonds(int nsb
, t_specbond sb
[])
125 for (i
= 0; (i
< nsb
); i
++)
131 sfree(sb
[i
].newres1
);
132 sfree(sb
[i
].newres2
);
136 static gmx_bool
is_special(int nsb
, t_specbond sb
[], char *res
, char *atom
)
140 for (i
= 0; (i
< nsb
); i
++)
142 if (((strncmp(sb
[i
].res1
, res
, 3) == 0) &&
143 (gmx_strcasecmp(sb
[i
].atom1
, atom
) == 0)) ||
144 ((strncmp(sb
[i
].res2
, res
, 3) == 0) &&
145 (gmx_strcasecmp(sb
[i
].atom2
, atom
) == 0)))
153 static gmx_bool
is_bond(int nsb
, t_specbond sb
[], t_atoms
*pdba
, int a1
, int a2
,
154 real d
, int *index_sb
, gmx_bool
*bSwap
)
157 char *at1
, *at2
, *res1
, *res2
;
159 at1
= *pdba
->atomname
[a1
];
160 at2
= *pdba
->atomname
[a2
];
161 res1
= *pdba
->resinfo
[pdba
->atom
[a1
].resind
].name
;
162 res2
= *pdba
->resinfo
[pdba
->atom
[a2
].resind
].name
;
166 fprintf(stderr
, "Checking %s-%d %s-%d and %s-%d %s-%d: %g ",
167 res1
, pdba
->resinfo
[pdba
->atom
[a1
].resind
].nr
, at1
, a1
+1,
168 res2
, pdba
->resinfo
[pdba
->atom
[a2
].resind
].nr
, at2
, a2
+1, d
);
171 for (i
= 0; (i
< nsb
); i
++)
174 if (((strncmp(sb
[i
].res1
, res1
, 3) == 0) &&
175 (gmx_strcasecmp(sb
[i
].atom1
, at1
) == 0) &&
176 (strncmp(sb
[i
].res2
, res2
, 3) == 0) &&
177 (gmx_strcasecmp(sb
[i
].atom2
, at2
) == 0)))
180 if ((0.9*sb
[i
].length
< d
) && (1.1*sb
[i
].length
> d
))
184 fprintf(stderr
, "%g\n", sb
[i
].length
);
189 if (((strncmp(sb
[i
].res1
, res2
, 3) == 0) &&
190 (gmx_strcasecmp(sb
[i
].atom1
, at2
) == 0) &&
191 (strncmp(sb
[i
].res2
, res1
, 3) == 0) &&
192 (gmx_strcasecmp(sb
[i
].atom2
, at1
) == 0)))
195 if ((0.9*sb
[i
].length
< d
) && (1.1*sb
[i
].length
> d
))
199 fprintf(stderr
, "%g\n", sb
[i
].length
);
207 fprintf(stderr
, "\n");
212 static void rename_1res(t_atoms
*pdba
, int resind
, char *newres
, gmx_bool bVerbose
)
216 printf("Using rtp entry %s for %s %d\n",
218 *pdba
->resinfo
[resind
].name
,
219 pdba
->resinfo
[resind
].nr
);
221 /* this used to free *resname, which messes up the symtab! */
222 snew(pdba
->resinfo
[resind
].rtp
, 1);
223 *pdba
->resinfo
[resind
].rtp
= gmx_strdup(newres
);
226 int mk_specbonds(t_atoms
*pdba
, rvec x
[], gmx_bool bInteractive
,
227 t_ssbond
**specbonds
, gmx_bool bVerbose
)
229 t_specbond
*sb
= nullptr;
230 t_ssbond
*bonds
= nullptr;
234 gmx_bool bDoit
, bSwap
;
236 int ai
, aj
, index_sb
;
241 sb
= get_specbonds(&nsb
);
245 snew(specp
, pdba
->nr
);
249 for (i
= 0; (i
< pdba
->nr
); i
++)
251 /* Check if this atom is special and if it is not a double atom
252 * in the input that still needs to be removed.
254 if (is_special(nsb
, sb
, *pdba
->resinfo
[pdba
->atom
[i
].resind
].name
,
255 *pdba
->atomname
[i
]) &&
257 pdba
->atom
[sgp
[nspec
-1]].resind
== pdba
->atom
[i
].resind
&&
258 gmx_strcasecmp(*pdba
->atomname
[sgp
[nspec
-1]],
259 *pdba
->atomname
[i
]) == 0))
261 specp
[nspec
] = pdba
->atom
[i
].resind
;
266 /* distance matrix d[nspec][nspec] */
268 for (i
= 0; (i
< nspec
); i
++)
273 for (i
= 0; (i
< nspec
); i
++)
276 for (j
= 0; (j
< nspec
); j
++)
279 d
[i
][j
] = std::sqrt(distance2(x
[ai
], x
[aj
]));
285 fprintf(stderr
, "Special Atom Distance matrix:\n");
286 for (b
= 0; (b
< nspec
); b
+= MAXCOL
)
288 /* print resname/number column headings */
289 fprintf(stderr
, "%8s%8s", "", "");
290 e
= std::min(b
+MAXCOL
, nspec
-1);
291 for (i
= b
; (i
< e
); i
++)
293 sprintf(buf
, "%s%d", *pdba
->resinfo
[pdba
->atom
[sgp
[i
]].resind
].name
,
294 pdba
->resinfo
[specp
[i
]].nr
);
295 fprintf(stderr
, "%8s", buf
);
297 fprintf(stderr
, "\n");
298 /* print atomname/number column headings */
299 fprintf(stderr
, "%8s%8s", "", "");
300 e
= std::min(b
+MAXCOL
, nspec
-1);
301 for (i
= b
; (i
< e
); i
++)
303 sprintf(buf
, "%s%d", *pdba
->atomname
[sgp
[i
]], sgp
[i
]+1);
304 fprintf(stderr
, "%8s", buf
);
306 fprintf(stderr
, "\n");
308 e
= std::min(b
+MAXCOL
, nspec
);
309 for (i
= b
+1; (i
< nspec
); i
++)
311 sprintf(buf
, "%s%d", *pdba
->resinfo
[pdba
->atom
[sgp
[i
]].resind
].name
,
312 pdba
->resinfo
[specp
[i
]].nr
);
313 fprintf(stderr
, "%8s", buf
);
314 sprintf(buf
, "%s%d", *pdba
->atomname
[sgp
[i
]], sgp
[i
]+1);
315 fprintf(stderr
, "%8s", buf
);
317 for (j
= b
; (j
< e2
); j
++)
319 fprintf(stderr
, " %7.3f", d
[i
][j
]);
321 fprintf(stderr
, "\n");
328 for (i
= 0; (i
< nspec
); i
++)
331 for (j
= i
+1; (j
< nspec
); j
++)
334 /* Ensure creation of at most nspec special bonds to avoid overflowing bonds[] */
335 if (nbonds
< nspec
&& is_bond(nsb
, sb
, pdba
, ai
, aj
, d
[i
][j
], &index_sb
, &bSwap
))
337 fprintf(stderr
, "%s %s-%d %s-%d and %s-%d %s-%d%s",
338 bInteractive
? "Link" : "Linking",
339 *pdba
->resinfo
[pdba
->atom
[ai
].resind
].name
,
340 pdba
->resinfo
[specp
[i
]].nr
,
341 *pdba
->atomname
[ai
], ai
+1,
342 *pdba
->resinfo
[pdba
->atom
[aj
].resind
].name
,
343 pdba
->resinfo
[specp
[j
]].nr
,
344 *pdba
->atomname
[aj
], aj
+1,
345 bInteractive
? " (y/n) ?" : "...\n");
346 bDoit
= bInteractive
? yesno() : TRUE
;
350 /* Store the residue numbers in the bonds array */
351 bonds
[nbonds
].res1
= specp
[i
];
352 bonds
[nbonds
].res2
= specp
[j
];
353 bonds
[nbonds
].a1
= gmx_strdup(*pdba
->atomname
[ai
]);
354 bonds
[nbonds
].a2
= gmx_strdup(*pdba
->atomname
[aj
]);
355 /* rename residues */
358 rename_1res(pdba
, specp
[i
], sb
[index_sb
].newres2
, bVerbose
);
359 rename_1res(pdba
, specp
[j
], sb
[index_sb
].newres1
, bVerbose
);
363 rename_1res(pdba
, specp
[i
], sb
[index_sb
].newres1
, bVerbose
);
364 rename_1res(pdba
, specp
[j
], sb
[index_sb
].newres2
, bVerbose
);
372 for (i
= 0; (i
< nspec
); i
++)
380 done_specbonds(nsb
, sb
);