?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixed printf format error
[gromacs.git] / src / kernel / resall.c
blobec3f932e81e37f9763a4426f89af0a48b6b94bbe
1 /*
2 * $Id$
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 3.2.0
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
15
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
20 *
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
27 *
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
30 *
31 * For more info, check our website at http://www.gromacs.org
32 *
33 * And Hey:
34 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
35 */
36 #include "sysstuff.h"
37 #include <ctype.h>
38 #include "assert.h"
39 #include "string2.h"
40 #include "strdb.h"
41 #include "futil.h"
42 #include "smalloc.h"
43 #include "fatal.h"
44 #include "symtab.h"
45 #include "macros.h"
46 #include "resall.h"
47 #include "pgutil.h"
48
49 t_atomtype *read_atype(char *adb,t_symtab *tab)
50 {
51 #define MAXAT 5000
52 FILE *in;
53 char aadb[STRLEN];
54 char buf[STRLEN],name[STRLEN];
55 double m;
56 int nratt;
57 t_atomtype *at;
58
59 sprintf(aadb,"%s.atp",adb);
60 in=libopen(aadb);
61
62 snew(at,1);
63 snew(at->atom,MAXAT);
64 snew(at->atomname,MAXAT);
65
66 for(nratt=0; ; nratt++) {
67 if (nratt >= MAXAT)
68 fatal_error(0,"nratt >= MAXAT(%d). Increase the latter",MAXAT);
69 if (feof(in))
70 break;
71
72 /* Skip blank or comment-only lines */
73 do {
74 fgets2(buf,STRLEN,in);
75 if(buf) {
76 strip_comment(buf);
77 trim(buf);
78 }
79 } while (buf && strlen(buf)==0);
80
81 if(buf==NULL)
82 break;
83
84 if (sscanf(buf,"%s%lf",name,&m) != 2)
85 break;
86 set_at(&(at->atom[nratt]),m,0.0,nratt,0);
87 at->atomname[nratt]=put_symtab(tab,name);
88 fprintf(stderr,"\rAtomtype %d",nratt+1);
89 }
90 fclose(in);
91 fprintf(stderr,"\n");
92 at->nr=nratt;
93
94 return at;
95 }
96
97 static void print_resatoms(FILE *out,t_atomtype *atype,t_restp *rtp)
98 {
99 int j,tp;
100
101 /* fprintf(out,"%5s\n",rtp->resname);
102 fprintf(out,"%5d\n",rtp->natom); */
103 fprintf(out,"[ %s ]\n",rtp->resname);
104 fprintf(out," [ atoms ]\n");
105
106 for(j=0; (j<rtp->natom); j++) {
107 tp=rtp->atom[j].type;
108 assert (tp >= 0);
109 assert (tp < atype->nr);
110 fprintf(out,"%6s%6s%8.3f%6d\n",
111 *(rtp->atomname[j]),*(atype->atomname[tp]),
112 rtp->atom[j].q,rtp->cgnr[j]);
113 }
114 }
115
116 static bool read_atoms(FILE *in,char *line,
117 t_restp *r0,t_symtab *tab,t_atomtype *atype)
118 {
119 int i,j,cg,maxentries;
120 char buf[256],buf1[256];
121 double q;
122
123 /* Read Atoms */
124 maxentries=0;
125 r0->atom= NULL;
126 r0->atomname= NULL;
127 r0->cgnr= NULL;
128 i=0;
129 while (get_a_line(in,line,STRLEN) && (strchr(line,'[')==NULL)) {
130 if (sscanf(line,"%s%s%lf%d",buf,buf1,&q,&cg) != 4)
131 return FALSE;
132 if (i>=maxentries) {
133 maxentries+=100;
134 srenew(r0->atom, maxentries);
135 srenew(r0->atomname, maxentries);
136 srenew(r0->cgnr, maxentries);
137 }
138 r0->atomname[i]=put_symtab(tab,buf);
139 r0->atom[i].q=q;
140 r0->cgnr[i]=cg;
141 for(j=0; (j<atype->nr); j++)
142 if (strcasecmp(buf1,*(atype->atomname[j])) == 0)
143 break;
144 if (j == atype->nr)
145 fatal_error(0,"Atom type %s (residue %s) not found in atomtype "
146 "database",buf1,r0->resname);
147 r0->atom[i].type=j;
148 r0->atom[i].m=atype->atom[j].m;
149 i++;
150 }
151 r0->natom=i;
152 srenew(r0->atom,i);
153 srenew(r0->atomname,i);
154 srenew(r0->cgnr,i);
155
156 return TRUE;
157 }
158
159 bool read_bondeds(int bt, FILE *in, char *line, t_restp *rtp)
160 {
161 char str[STRLEN];
162 int j,n,ni,maxrb;
163
164 maxrb = rtp->rb[bt].nb;
165 while (get_a_line(in,line,STRLEN) && (strchr(line,'[')==NULL)) {
166 if ( rtp->rb[bt].nb >= maxrb ) {
167 maxrb+=100;
168 srenew(rtp->rb[bt].b,maxrb);
169 }
170 n=0;
171 for(j=0; j < btsNiatoms[bt]; j++) {
172 if ( sscanf(line+n,"%s%n",str,&ni)==1 )
173 rtp->rb[bt].b[rtp->rb[bt].nb].a[j]=strdup(str);
174 else
175 return FALSE;
176 n+=ni;
177 }
178 for( ; j < MAXATOMLIST; j++)
179 rtp->rb[bt].b[rtp->rb[bt].nb].a[j]=NULL;
180 while (isspace(line[n]))
181 n++;
182 rtrim(line+n);
183 rtp->rb[bt].b[rtp->rb[bt].nb].s=strdup(line+n);
184 rtp->rb[bt].nb++;
185 }
186 /* give back unused memory */
187 srenew(rtp->rb[bt].b,rtp->rb[bt].nb);
188
189 return TRUE;
190 }
191
192 static void print_resbondeds(FILE *out, int bt, t_restp *rtp)
193 {
194 int i,j;
195
196 if (rtp->rb[bt].nb) {
197 fprintf(out," [ %s ]\n",btsNames[bt]);
198
199 for(i=0; i < rtp->rb[bt].nb; i++) {
200 for(j=0; j<btsNiatoms[bt]; j++)
201 fprintf(out,"%6s ",rtp->rb[bt].b[i].a[j]);
202 if (rtp->rb[bt].b[i].s[0])
203 fprintf(out," %s",rtp->rb[bt].b[i].s);
204 fprintf(out,"\n");
205 }
206 }
207 }
208
209 static void check_rtp(int nrtp,t_restp rtp[],char *libfn)
210 {
211 int i;
212
213 /* check for double entries, assuming list is already sorted */
214 for(i=1; (i<nrtp); i++) {
215 if (strcasecmp(rtp[i-1].resname,rtp[i].resname) == 0)
216 fprintf(stderr,"WARNING double entry %s in file %s\n",
217 rtp[i].resname,libfn);
218 }
219 }
220
221 static int comprtp(const void *a,const void *b)
222 {
223 t_restp *ra,*rb;
224
225 ra=(t_restp *)a;
226 rb=(t_restp *)b;
227
228 return strcasecmp(ra->resname,rb->resname);
229 }
230
231 int get_bt(char* header)
232 {
233 int i;
234
235 for(i=0; i<ebtsNR; i++)
236 if ( strcasecmp(btsNames[i],header)==0 )
237 return i;
238 return NOTSET;
239 }
240
241 void clear_t_restp(t_restp *rrtp)
242 {
243 memset((void *)rrtp, 0, sizeof(t_restp));
244 }
245
246 int read_resall(char *ff, int bts[], t_restp **rtp,
247 t_atomtype *atype, t_symtab *tab, bool *bAlldih, int *nrexcl,
248 bool *HH14, bool *bRemoveDih)
249 {
250 FILE *in;
251 char rrdb[STRLEN],line[STRLEN],header[STRLEN];
252 int i,nrtp,maxrtp,bt,nparam;
253 t_restp *rrtp;
254 bool bNextResidue,bError;
255
256 sprintf(rrdb,"%s.rtp",ff);
257 in=libopen(rrdb);
258 rrtp=NULL;
259 if (debug) {
260 fprintf(debug,"%9s %5s", "Residue", "atoms");
261 for(i=0; i<ebtsNR; i++)
262 fprintf(debug," %10s",btsNames[i]);
263 fprintf(debug,"\n");
264 }
265
266 /* these bonded parameters will overwritten be when *
267 * there is a [ bondedtypes ] entry in the .rtp file */
268 bts[0] = 1; /* normal bonds */
269 bts[1] = 1; /* normal angles */
270 bts[2] = 1; /* normal dihedrals */
271 bts[3] = 2; /* normal impropers */
272
273
274 /* Column 5 & 6 aren't really bonded types, but we include
275 * them here to avoid introducing a new section:
276 * Column 5: 1 means generate all dihedrals, 0 not.
277 * Column 6: Number of bonded neighbors to exclude.
278 * Coulmn 7: Generate 1,4 interactions between pairs of hydrogens
279 * Column 8: Remove impropers over the same bond as a proper dihedral
280 */
281
282 nrtp=0;
283 maxrtp=0;
284 get_a_line(in,line,STRLEN);
285 if (!get_header(line,header))
286 fatal_error(0,"in .rtp file at line:\n%s\n",line);
287 if (strncasecmp("bondedtypes",header,5)==0) {
288 get_a_line(in,line,STRLEN);
289 if ((nparam=sscanf(line,"%d %d %d %d %d %d %d %d",&bts[0],&bts[1],&bts[2],&bts[3],bAlldih,nrexcl,HH14,bRemoveDih)) < 4 )
290 fatal_error(0,"need at least 4 (up to 8) parameters in .rtp file at line:\n%s\n",line);
291 get_a_line(in,line,STRLEN);
292 if(nparam<5) {
293 fprintf(stderr,"Using default: not generating all possible dihedrals\n");
294 *bAlldih=FALSE;
295 }
296 if(nparam<6) {
297 fprintf(stderr,"Using default: excluding 3 bonded neighbors\n");
298 *nrexcl=3;
299 }
300 if(nparam<7) {
301 fprintf(stderr,"Using default: generating 1,4 H--H interactions\n");
302 *HH14=TRUE;
303 }
304 if(nparam<8) {
305 fprintf(stderr,"Using default: removing impropers on same bond as a proper\n");
306 *bRemoveDih=TRUE;
307 }
308 } else {
309 fprintf(stderr,
310 "Reading .rtp file without '[ bondedtypes ]' directive,\n"
311 "Will proceed as if the entry\n"
312 "\n"
313 "\n[ bondedtypes ]"
314 "\n; bonds angles dihedrals impropers all_dihedrals nr_exclusions HH14 remove_dih"
315 "\n %3d %3d %3d %3d %3d %3d %3d %3d"
316 "\n"
317 "was present at the beginning of %s",
318 bts[0],bts[1],bts[2],bts[3], (*bAlldih) ? 1 : 0,*nrexcl,*HH14,*bRemoveDih,rrdb);
319 }
320 nrtp=0;
321 while (!feof(in)) {
322 if (nrtp >= maxrtp) {
323 maxrtp+=100;
324 srenew(rrtp,maxrtp);
325 }
326 clear_t_restp(&rrtp[nrtp]);
327 if (!get_header(line,header))
328 fatal_error(0,"in .rtp file at line:\n%s\n",line);
329 rrtp[nrtp].resname=strdup(header);
330
331 get_a_line(in,line,STRLEN);
332 bError=FALSE;
333 bNextResidue=FALSE;
334 while (!bNextResidue) {
335 if (!get_header(line,header))
336 if (feof(in))
337 bNextResidue=TRUE;
338 else
339 bError=TRUE;
340 else {
341 bt=get_bt(header);
342 if ( bt!=NOTSET )
343 bError=!read_bondeds(bt,in,line,&rrtp[nrtp]);
344 else
345 if (strncasecmp("atoms",header,5)==0)
346 bError=!read_atoms(in,line,&(rrtp[nrtp]),tab,atype);
347 else {
348 if (!feof(in) && !get_header(line,header))
349 bError=TRUE;
350 bNextResidue=TRUE;
351 }
352 }
353 if (bError)
354 fatal_error(0,"in .rtp file in residue %s at line:\n%s\n",
355 rrtp[nrtp].resname,line);
356 }
357 if (rrtp[nrtp].natom == 0)
358 fatal_error(0,"No atoms found in .rtp file in residue %s\n",
359 rrtp[nrtp].resname);
360 if (debug) {
361 fprintf(debug,"%3d %5s %5d",
362 nrtp+1,rrtp[nrtp].resname,rrtp[nrtp].natom);
363 for(i=0; i<ebtsNR; i++)
364 fprintf(debug," %10d",rrtp[nrtp].rb[i].nb);
365 fprintf(debug,"\n");
366 }
367 nrtp++;
368 fprintf(stderr,"\rResidue %d",nrtp);
369 }
370 fclose(in);
371 /* give back unused memory */
372 srenew(rrtp,nrtp);
373
374 fprintf(stderr,"\nSorting it all out...\n");
375 qsort(rrtp,nrtp,(size_t)sizeof(rrtp[0]),comprtp);
376
377 check_rtp(nrtp,rrtp,rrdb);
378
379 *rtp = rrtp;
380
381 return nrtp;
382 }
383
384 void print_resall(FILE *out, int bts[], int nrtp, t_restp rtp[],
385 t_atomtype *atype, bool bAlldih, int nrexcl,
386 bool HH14, bool bRemoveDih)
387 {
388 int i,bt;
389
390 /* print all the ebtsNR type numbers */
391 fprintf(out,"[ bondedtypes ]\n");
392 fprintf(out,"; bonds angles dihedrals impropers all_dihedrals nr_exclusions HH14 remove_dih\n");
393 fprintf(out," %5d %6d %9d %9d %14d %14d %14d %14d\n\n",bts[0],bts[1],bts[2],bts[3],bAlldih,nrexcl,HH14,bRemoveDih);
394
395 for(i=0; i<nrtp; i++) {
396 if (rtp[i].natom > 0) {
397 print_resatoms(out,atype,&rtp[i]);
398 for(bt=0; bt<ebtsNR; bt++)
399 print_resbondeds(out,bt,&rtp[i]);
400 }
401 }
402 }
403
404 /************************************************************
405 *
406 * SEARCH ROUTINES
407 *
408 ***********************************************************/
409 int neq_str(char *a1,char *a2)
410 {
411 int j,l;
412
413 l=min((int)strlen(a1),(int)strlen(a2));
414 j=0;
415 while ( (j<l) && (toupper(a1[j]) == toupper(a2[j])) )
416 j++;
417
418 return j;
419 }
420
421 t_restp *search_rtp(char *key,int nrtp,t_restp rtp[])
422 {
423 int i,n,best,besti;
424
425 besti=-1;
426 best=1;
427 for(i=0; (i<nrtp); i++) {
428 n=neq_str(key,rtp[i].resname);
429 if (n > best) {
430 besti=i;
431 best=n;
432 }
433 }
434 if (besti == -1)
435 fatal_error(0,"Residue '%s' not found in residue topology database\n",key);
436 if (strlen(rtp[besti].resname) != strlen(key))
437 fprintf(stderr,"Warning: '%s' not found in residue topology database, "
438 "trying to use '%s'\n", key, rtp[besti].resname);
439
440 return &rtp[besti];
441 }
The ultimate molecular dynamics simulation package