3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
55 c=toupper(fgetc(stdin));
56 } while ((c != 'Y') && (c != 'N'));
61 t_specbond *get_specbonds(int *nspecbond)
63 const char *sbfile="specbond.dat";
66 char r1buf[32],r2buf[32],a1buf[32],a2buf[32],nr1buf[32],nr2buf[32];
72 nlines = get_lines(sbfile,&lines);
77 for(i=0; (i<nlines); i++) {
78 if (sscanf(lines[i],"%s%s%d%s%s%d%lf%s%s",
79 r1buf,a1buf,&nb1,r2buf,a2buf,&nb2,&length,nr1buf,nr2buf) != 9)
80 fprintf(stderr,"Invalid line '%s' in %s\n",lines[i],sbfile);
82 sb[n].res1 = strdup(r1buf);
83 sb[n].res2 = strdup(r2buf);
84 sb[n].newres1= strdup(nr1buf);
85 sb[n].newres2= strdup(nr2buf);
86 sb[n].atom1 = strdup(a1buf);
87 sb[n].atom2 = strdup(a2buf);
90 sb[n].length = length;
97 fprintf(stderr,"%d out of %d lines of %s converted successfully\n",
105 void done_specbonds(int nsb,t_specbond sb[])
109 for(i=0; (i<nsb); i++) {
114 sfree(sb[i].newres1);
115 sfree(sb[i].newres2);
119 static bool is_special(int nsb,t_specbond sb[],char *res,char *atom)
123 for(i=0; (i<nsb); i++) {
124 if (((strncmp(sb[i].res1,res,3) == 0) &&
125 (gmx_strcasecmp(sb[i].atom1,atom) == 0)) ||
126 ((strncmp(sb[i].res2,res,3) == 0) &&
127 (gmx_strcasecmp(sb[i].atom2,atom) == 0)))
133 static bool is_bond(int nsb,t_specbond sb[],t_atoms *pdba,int a1,int a2,
134 real d,int *index_sb,bool *bSwap)
137 char *at1,*at2,*res1,*res2;
139 at1=*pdba->atomname[a1];
140 at2=*pdba->atomname[a2];
141 res1=*pdba->resinfo[pdba->atom[a1].resind].name;
142 res2=*pdba->resinfo[pdba->atom[a2].resind].name;
145 fprintf(stderr,"Checking %s-%d %s-%d and %s-%d %s-%d: %g ",
146 res1, pdba->resinfo[pdba->atom[a1].resind].nr, at1, a1+1,
147 res2, pdba->resinfo[pdba->atom[a2].resind].nr, at2, a2+1, d);
149 for(i=0; (i<nsb); i++) {
151 if (((strncmp(sb[i].res1,res1,3) == 0) &&
152 (gmx_strcasecmp(sb[i].atom1,at1) == 0) &&
153 (strncmp(sb[i].res2,res2,3) == 0) &&
154 (gmx_strcasecmp(sb[i].atom2,at2) == 0))) {
156 if ((0.9*sb[i].length < d) && (1.1*sb[i].length > d)) {
157 if (debug) fprintf(stderr,"%g\n", sb[i].length);
161 if (((strncmp(sb[i].res1,res2,3) == 0) &&
162 (gmx_strcasecmp(sb[i].atom1,at2) == 0) &&
163 (strncmp(sb[i].res2,res1,3) == 0) &&
164 (gmx_strcasecmp(sb[i].atom2,at1) == 0))) {
166 if ((0.9*sb[i].length < d) && (1.1*sb[i].length > d)) {
167 if (debug) fprintf(stderr,"%g\n", sb[i].length);
172 if (debug) fprintf(stderr,"\n");
176 static void rename_1res(t_atoms *pdba,int resind,char *newres,bool bVerbose)
179 printf("Using rtp entry %s for %s %d\n",
181 *pdba->resinfo[resind].name,
182 pdba->resinfo[resind].nr);
184 /* this used to free *resname, which fucks up the symtab! */
185 snew(pdba->resinfo[resind].rtp,1);
186 *pdba->resinfo[resind].rtp = strdup(newres);
189 int mk_specbonds(t_atoms *pdba,rvec x[],bool bInteractive,
190 t_ssbond **specbonds,bool bVerbose)
193 t_ssbond *bonds=NULL;
204 sb = get_specbonds(&nsb);
207 snew(specp,pdba->nr);
211 for(i=0;(i<pdba->nr);i++) {
212 /* Check if this atom is special and if it is not a double atom
213 * in the input that still needs to be removed.
215 if (is_special(nsb,sb,*pdba->resinfo[pdba->atom[i].resind].name,
216 *pdba->atomname[i]) &&
218 pdba->atom[sgp[nspec-1]].resind == pdba->atom[i].resind &&
219 gmx_strcasecmp(*pdba->atomname[sgp[nspec-1]],
220 *pdba->atomname[i]) == 0)) {
221 specp[nspec] = pdba->atom[i].resind;
226 /* distance matrix d[nspec][nspec] */
228 for(i=0; (i<nspec); i++)
231 for(i=0; (i<nspec); i++) {
233 for(j=0; (j<nspec); j++) {
235 d[i][j]=sqrt(distance2(x[ai],x[aj]));
240 fprintf(stderr,"Special Atom Distance matrix:\n");
241 for(b=0; (b<nspec); b+=MAXCOL) {
242 /* print resname/number column headings */
243 fprintf(stderr,"%8s%8s","","");
244 e=min(b+MAXCOL, nspec-1);
245 for(i=b; (i<e); i++) {
246 sprintf(buf,"%s%d",*pdba->resinfo[pdba->atom[sgp[i]].resind].name,
247 pdba->resinfo[specp[i]].nr);
248 fprintf(stderr,"%8s",buf);
250 fprintf(stderr,"\n");
251 /* print atomname/number column headings */
252 fprintf(stderr,"%8s%8s","","");
253 e=min(b+MAXCOL, nspec-1);
254 for(i=b; (i<e); i++) {
255 sprintf(buf,"%s%d",*pdba->atomname[sgp[i]], sgp[i]+1);
256 fprintf(stderr,"%8s",buf);
258 fprintf(stderr,"\n");
260 e=min(b+MAXCOL, nspec);
261 for(i=b+1; (i<nspec); i++) {
262 sprintf(buf,"%s%d",*pdba->resinfo[pdba->atom[sgp[i]].resind].name,
263 pdba->resinfo[specp[i]].nr);
264 fprintf(stderr,"%8s",buf);
265 sprintf(buf,"%s%d", *pdba->atomname[sgp[i]], sgp[i]+1);
266 fprintf(stderr,"%8s",buf);
268 for(j=b; (j<e2); j++)
269 fprintf(stderr," %7.3f",d[i][j]);
270 fprintf(stderr,"\n");
277 for(i=0; (i<nspec); i++) {
279 for(j=i+1; (j<nspec); j++) {
281 if (is_bond(nsb,sb,pdba,ai,aj,d[i][j],&index_sb,&bSwap)) {
282 fprintf(stderr,"%s %s-%d %s-%d and %s-%d %s-%d%s",
283 bInteractive ? "Link" : "Linking",
284 *pdba->resinfo[pdba->atom[ai].resind].name,
285 pdba->resinfo[specp[i]].nr,
286 *pdba->atomname[ai], ai+1,
287 *pdba->resinfo[pdba->atom[aj].resind].name,
288 pdba->resinfo[specp[j]].nr,
289 *pdba->atomname[aj], aj+1,
290 bInteractive ? " (y/n) ?" : "...\n");
291 bDoit=bInteractive ? yesno() : TRUE;
294 /* Store the residue numbers in the bonds array */
295 bonds[nbonds].res1 = specp[i];
296 bonds[nbonds].res2 = specp[j];
297 bonds[nbonds].a1 = strdup(*pdba->atomname[ai]);
298 bonds[nbonds].a2 = strdup(*pdba->atomname[aj]);
299 /* rename residues */
301 rename_1res(pdba,specp[i],sb[index_sb].newres2,bVerbose);
302 rename_1res(pdba,specp[j],sb[index_sb].newres1,bVerbose);
305 rename_1res(pdba,specp[i],sb[index_sb].newres1,bVerbose);
306 rename_1res(pdba,specp[j],sb[index_sb].newres2,bVerbose);
314 for(i=0; (i<nspec); i++)
320 done_specbonds(nsb,sb);