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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
58 c=toupper(fgetc(stdin));
59 } while ((c != 'Y') && (c != 'N'));
64 t_specbond *get_specbonds(int *nspecbond)
66 const char *sbfile="specbond.dat";
69 char r1buf[32],r2buf[32],a1buf[32],a2buf[32],nr1buf[32],nr2buf[32];
75 nlines = get_lines(sbfile,&lines);
80 for(i=0; (i<nlines); i++) {
81 if (sscanf(lines[i],"%s%s%d%s%s%d%lf%s%s",
82 r1buf,a1buf,&nb1,r2buf,a2buf,&nb2,&length,nr1buf,nr2buf) != 9)
83 fprintf(stderr,"Invalid line '%s' in %s\n",lines[i],sbfile);
85 sb[n].res1 = strdup(r1buf);
86 sb[n].res2 = strdup(r2buf);
87 sb[n].newres1= strdup(nr1buf);
88 sb[n].newres2= strdup(nr2buf);
89 sb[n].atom1 = strdup(a1buf);
90 sb[n].atom2 = strdup(a2buf);
93 sb[n].length = length;
100 fprintf(stderr,"%d out of %d lines of %s converted successfully\n",
108 void done_specbonds(int nsb,t_specbond sb[])
112 for(i=0; (i<nsb); i++) {
117 sfree(sb[i].newres1);
118 sfree(sb[i].newres2);
122 static gmx_bool is_special(int nsb,t_specbond sb[],char *res,char *atom)
126 for(i=0; (i<nsb); i++) {
127 if (((strncmp(sb[i].res1,res,3) == 0) &&
128 (gmx_strcasecmp(sb[i].atom1,atom) == 0)) ||
129 ((strncmp(sb[i].res2,res,3) == 0) &&
130 (gmx_strcasecmp(sb[i].atom2,atom) == 0)))
136 static gmx_bool is_bond(int nsb,t_specbond sb[],t_atoms *pdba,int a1,int a2,
137 real d,int *index_sb,gmx_bool *bSwap)
140 char *at1,*at2,*res1,*res2;
142 at1=*pdba->atomname[a1];
143 at2=*pdba->atomname[a2];
144 res1=*pdba->resinfo[pdba->atom[a1].resind].name;
145 res2=*pdba->resinfo[pdba->atom[a2].resind].name;
148 fprintf(stderr,"Checking %s-%d %s-%d and %s-%d %s-%d: %g ",
149 res1, pdba->resinfo[pdba->atom[a1].resind].nr, at1, a1+1,
150 res2, pdba->resinfo[pdba->atom[a2].resind].nr, at2, a2+1, d);
152 for(i=0; (i<nsb); i++) {
154 if (((strncmp(sb[i].res1,res1,3) == 0) &&
155 (gmx_strcasecmp(sb[i].atom1,at1) == 0) &&
156 (strncmp(sb[i].res2,res2,3) == 0) &&
157 (gmx_strcasecmp(sb[i].atom2,at2) == 0))) {
159 if ((0.9*sb[i].length < d) && (1.1*sb[i].length > d)) {
160 if (debug) fprintf(stderr,"%g\n", sb[i].length);
164 if (((strncmp(sb[i].res1,res2,3) == 0) &&
165 (gmx_strcasecmp(sb[i].atom1,at2) == 0) &&
166 (strncmp(sb[i].res2,res1,3) == 0) &&
167 (gmx_strcasecmp(sb[i].atom2,at1) == 0))) {
169 if ((0.9*sb[i].length < d) && (1.1*sb[i].length > d)) {
170 if (debug) fprintf(stderr,"%g\n", sb[i].length);
175 if (debug) fprintf(stderr,"\n");
179 static void rename_1res(t_atoms *pdba,int resind,char *newres,gmx_bool bVerbose)
182 printf("Using rtp entry %s for %s %d\n",
184 *pdba->resinfo[resind].name,
185 pdba->resinfo[resind].nr);
187 /* this used to free *resname, which fucks up the symtab! */
188 snew(pdba->resinfo[resind].rtp,1);
189 *pdba->resinfo[resind].rtp = strdup(newres);
192 int mk_specbonds(t_atoms *pdba,rvec x[],gmx_bool bInteractive,
193 t_ssbond **specbonds,gmx_bool bVerbose)
196 t_ssbond *bonds=NULL;
200 gmx_bool bDoit,bSwap;
207 sb = get_specbonds(&nsb);
210 snew(specp,pdba->nr);
214 for(i=0;(i<pdba->nr);i++) {
215 /* Check if this atom is special and if it is not a double atom
216 * in the input that still needs to be removed.
218 if (is_special(nsb,sb,*pdba->resinfo[pdba->atom[i].resind].name,
219 *pdba->atomname[i]) &&
221 pdba->atom[sgp[nspec-1]].resind == pdba->atom[i].resind &&
222 gmx_strcasecmp(*pdba->atomname[sgp[nspec-1]],
223 *pdba->atomname[i]) == 0)) {
224 specp[nspec] = pdba->atom[i].resind;
229 /* distance matrix d[nspec][nspec] */
231 for(i=0; (i<nspec); i++)
234 for(i=0; (i<nspec); i++) {
236 for(j=0; (j<nspec); j++) {
238 d[i][j]=sqrt(distance2(x[ai],x[aj]));
243 fprintf(stderr,"Special Atom Distance matrix:\n");
244 for(b=0; (b<nspec); b+=MAXCOL) {
245 /* print resname/number column headings */
246 fprintf(stderr,"%8s%8s","","");
247 e=min(b+MAXCOL, nspec-1);
248 for(i=b; (i<e); i++) {
249 sprintf(buf,"%s%d",*pdba->resinfo[pdba->atom[sgp[i]].resind].name,
250 pdba->resinfo[specp[i]].nr);
251 fprintf(stderr,"%8s",buf);
253 fprintf(stderr,"\n");
254 /* print atomname/number column headings */
255 fprintf(stderr,"%8s%8s","","");
256 e=min(b+MAXCOL, nspec-1);
257 for(i=b; (i<e); i++) {
258 sprintf(buf,"%s%d",*pdba->atomname[sgp[i]], sgp[i]+1);
259 fprintf(stderr,"%8s",buf);
261 fprintf(stderr,"\n");
263 e=min(b+MAXCOL, nspec);
264 for(i=b+1; (i<nspec); i++) {
265 sprintf(buf,"%s%d",*pdba->resinfo[pdba->atom[sgp[i]].resind].name,
266 pdba->resinfo[specp[i]].nr);
267 fprintf(stderr,"%8s",buf);
268 sprintf(buf,"%s%d", *pdba->atomname[sgp[i]], sgp[i]+1);
269 fprintf(stderr,"%8s",buf);
271 for(j=b; (j<e2); j++)
272 fprintf(stderr," %7.3f",d[i][j]);
273 fprintf(stderr,"\n");
280 for(i=0; (i<nspec); i++) {
282 for(j=i+1; (j<nspec); j++) {
284 if (is_bond(nsb,sb,pdba,ai,aj,d[i][j],&index_sb,&bSwap)) {
285 fprintf(stderr,"%s %s-%d %s-%d and %s-%d %s-%d%s",
286 bInteractive ? "Link" : "Linking",
287 *pdba->resinfo[pdba->atom[ai].resind].name,
288 pdba->resinfo[specp[i]].nr,
289 *pdba->atomname[ai], ai+1,
290 *pdba->resinfo[pdba->atom[aj].resind].name,
291 pdba->resinfo[specp[j]].nr,
292 *pdba->atomname[aj], aj+1,
293 bInteractive ? " (y/n) ?" : "...\n");
294 bDoit=bInteractive ? yesno() : TRUE;
297 /* Store the residue numbers in the bonds array */
298 bonds[nbonds].res1 = specp[i];
299 bonds[nbonds].res2 = specp[j];
300 bonds[nbonds].a1 = strdup(*pdba->atomname[ai]);
301 bonds[nbonds].a2 = strdup(*pdba->atomname[aj]);
302 /* rename residues */
304 rename_1res(pdba,specp[i],sb[index_sb].newres2,bVerbose);
305 rename_1res(pdba,specp[j],sb[index_sb].newres1,bVerbose);
308 rename_1res(pdba,specp[i],sb[index_sb].newres1,bVerbose);
309 rename_1res(pdba,specp[j],sb[index_sb].newres2,bVerbose);
317 for(i=0; (i<nspec); i++)
323 done_specbonds(nsb,sb);