1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
31 * For more info, check our website at http://www.gromacs.org
34 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
46 #include "gmx_fatal.h"
51 #include "fflibutil.h"
53 gpp_atomtype_t
read_atype(const char *ffdir
,t_symtab
*tab
)
58 char buf
[STRLEN
],name
[STRLEN
];
65 nfile
= fflib_search_file_end(ffdir
,".atp",TRUE
,&file
);
70 for(f
=0; f
<nfile
; f
++)
72 in
= fflib_open(file
[f
]);
75 /* Skip blank or comment-only lines */
78 fgets2(buf
,STRLEN
,in
);
85 while (!feof(in
) && NULL
!=buf
&& strlen(buf
)==0);
87 if ((buf
!= NULL
) && (sscanf(buf
,"%s%lf",name
,&m
) == 2))
90 add_atomtype(at
,tab
,a
,name
,nb
, 0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 );
91 fprintf(stderr
,"\rAtomtype %d",nratt
+1);
103 static void print_resatoms(FILE *out
,gpp_atomtype_t atype
,t_restp
*rtp
)
108 /* fprintf(out,"%5s\n",rtp->resname);
109 fprintf(out,"%5d\n",rtp->natom); */
110 fprintf(out
,"[ %s ]\n",rtp
->resname
);
111 fprintf(out
," [ atoms ]\n");
113 for(j
=0; (j
<rtp
->natom
); j
++) {
114 tp
= rtp
->atom
[j
].type
;
115 tpnm
= get_atomtype_name(tp
,atype
);
117 gmx_fatal(FARGS
,"Incorrect atomtype (%d)",tp
);
118 fprintf(out
,"%6s %6s %8.3f %6d\n",
119 *(rtp
->atomname
[j
]),tpnm
,rtp
->atom
[j
].q
,rtp
->cgnr
[j
]);
123 static gmx_bool
read_atoms(FILE *in
,char *line
,
124 t_restp
*r0
,t_symtab
*tab
,gpp_atomtype_t atype
)
126 int i
,j
,cg
,maxentries
;
127 char buf
[256],buf1
[256];
136 while (get_a_line(in
,line
,STRLEN
) && (strchr(line
,'[')==NULL
)) {
137 if (sscanf(line
,"%s%s%lf%d",buf
,buf1
,&q
,&cg
) != 4)
141 srenew(r0
->atom
, maxentries
);
142 srenew(r0
->atomname
, maxentries
);
143 srenew(r0
->cgnr
, maxentries
);
145 r0
->atomname
[i
] = put_symtab(tab
,buf
);
148 j
= get_atomtype_type(buf1
,atype
);
150 gmx_fatal(FARGS
,"Atom type %s (residue %s) not found in atomtype "
151 "database",buf1
,r0
->resname
);
153 r0
->atom
[i
].m
=get_atomtype_massA(j
,atype
);
158 srenew(r0
->atomname
,i
);
164 gmx_bool
read_bondeds(int bt
, FILE *in
, char *line
, t_restp
*rtp
)
169 maxrb
= rtp
->rb
[bt
].nb
;
170 while (get_a_line(in
,line
,STRLEN
) && (strchr(line
,'[')==NULL
)) {
171 if ( rtp
->rb
[bt
].nb
>= maxrb
) {
173 srenew(rtp
->rb
[bt
].b
,maxrb
);
176 for(j
=0; j
< btsNiatoms
[bt
]; j
++) {
177 if ( sscanf(line
+n
,"%s%n",str
,&ni
)==1 )
178 rtp
->rb
[bt
].b
[rtp
->rb
[bt
].nb
].a
[j
]=strdup(str
);
183 for( ; j
< MAXATOMLIST
; j
++)
184 rtp
->rb
[bt
].b
[rtp
->rb
[bt
].nb
].a
[j
]=NULL
;
185 while (isspace(line
[n
]))
188 rtp
->rb
[bt
].b
[rtp
->rb
[bt
].nb
].s
=strdup(line
+n
);
191 /* give back unused memory */
192 srenew(rtp
->rb
[bt
].b
,rtp
->rb
[bt
].nb
);
197 static void print_resbondeds(FILE *out
, int bt
, t_restp
*rtp
)
201 if (rtp
->rb
[bt
].nb
) {
202 fprintf(out
," [ %s ]\n",btsNames
[bt
]);
204 for(i
=0; i
< rtp
->rb
[bt
].nb
; i
++) {
205 for(j
=0; j
<btsNiatoms
[bt
]; j
++)
206 fprintf(out
,"%6s ",rtp
->rb
[bt
].b
[i
].a
[j
]);
207 if (rtp
->rb
[bt
].b
[i
].s
[0])
208 fprintf(out
," %s",rtp
->rb
[bt
].b
[i
].s
);
214 static void check_rtp(int nrtp
,t_restp rtp
[],char *libfn
)
218 /* check for double entries, assuming list is already sorted */
219 for(i
=1; (i
<nrtp
); i
++) {
220 if (gmx_strcasecmp(rtp
[i
-1].resname
,rtp
[i
].resname
) == 0)
221 fprintf(stderr
,"WARNING double entry %s in file %s\n",
222 rtp
[i
].resname
,libfn
);
226 static int comprtp(const void *a
,const void *b
)
233 return gmx_strcasecmp(ra
->resname
,rb
->resname
);
236 int get_bt(char* header
)
240 for(i
=0; i
<ebtsNR
; i
++)
241 if ( gmx_strcasecmp(btsNames
[i
],header
)==0 )
246 void clear_t_restp(t_restp
*rrtp
)
248 memset((void *)rrtp
, 0, sizeof(t_restp
));
251 /* print all the ebtsNR type numbers */
252 void print_resall_header(FILE *out
, t_restp rtp
[])
254 fprintf(out
,"[ bondedtypes ]\n");
255 fprintf(out
,"; bonds angles dihedrals impropers all_dihedrals nr_exclusions HH14 remove_dih\n");
256 fprintf(out
," %5d %6d %9d %9d %14d %14d %14d %14d\n\n",
261 rtp
[0].bKeepAllGeneratedDihedrals
,
263 rtp
[0].bGenerateHH14Interactions
,
264 rtp
[0].bRemoveDihedralIfWithImproper
);
267 void print_resall(FILE *out
, int nrtp
, t_restp rtp
[],
268 gpp_atomtype_t atype
)
276 print_resall_header(out
, rtp
);
278 for(i
=0; i
<nrtp
; i
++) {
279 if (rtp
[i
].natom
> 0) {
280 print_resatoms(out
,atype
,&rtp
[i
]);
281 for(bt
=0; bt
<ebtsNR
; bt
++)
282 print_resbondeds(out
,bt
,&rtp
[i
]);
287 void read_resall(char *rrdb
, int *nrtpptr
, t_restp
**rtp
,
288 gpp_atomtype_t atype
, t_symtab
*tab
,
289 gmx_bool bAllowOverrideRTP
)
292 char filebase
[STRLEN
],*ptr
,line
[STRLEN
],header
[STRLEN
];
293 int i
,nrtp
,maxrtp
,bt
,nparam
;
295 t_restp
*rrtp
, *header_settings
;
296 gmx_bool bNextResidue
,bError
;
299 fflib_filename_base(rrdb
,filebase
,STRLEN
);
301 in
= fflib_open(rrdb
);
304 fprintf(debug
,"%9s %5s", "Residue", "atoms");
305 for(i
=0; i
<ebtsNR
; i
++)
306 fprintf(debug
," %10s",btsNames
[i
]);
309 snew(header_settings
, 1);
311 /* these bonded parameters will overwritten be when *
312 * there is a [ bondedtypes ] entry in the .rtp file */
313 header_settings
->rb
[ebtsBONDS
].type
= 1; /* normal bonds */
314 header_settings
->rb
[ebtsANGLES
].type
= 1; /* normal angles */
315 header_settings
->rb
[ebtsPDIHS
].type
= 1; /* normal dihedrals */
316 header_settings
->rb
[ebtsIDIHS
].type
= 2; /* normal impropers */
317 header_settings
->rb
[ebtsEXCLS
].type
= 1; /* normal exclusions */
318 header_settings
->rb
[ebtsCMAP
].type
= 1; /* normal cmap torsions */
320 header_settings
->bKeepAllGeneratedDihedrals
= FALSE
;
321 header_settings
->nrexcl
= 3;
322 header_settings
->bGenerateHH14Interactions
= TRUE
;
323 header_settings
->bRemoveDihedralIfWithImproper
= TRUE
;
325 /* Column 5 & 6 aren't really bonded types, but we include
326 * them here to avoid introducing a new section:
327 * Column 5 : This controls the generation of dihedrals from the bonding.
328 * All possible dihedrals are generated automatically. A value of
329 * 1 here means that all these are retained. A value of
330 * 0 here requires generated dihedrals be removed if
331 * * there are any dihedrals on the same central atoms
332 * specified in the residue topology, or
333 * * there are other identical generated dihedrals
334 * sharing the same central atoms, or
335 * * there are other generated dihedrals sharing the
336 * same central bond that have fewer hydrogen atoms
337 * Column 6: Number of bonded neighbors to exclude.
338 * Column 7: Generate 1,4 interactions between two hydrogen atoms
339 * Column 8: Remove proper dihedrals if centered on the same bond
340 * as an improper dihedral
342 get_a_line(in
,line
,STRLEN
);
343 if (!get_header(line
,header
))
344 gmx_fatal(FARGS
,"in .rtp file at line:\n%s\n",line
);
345 if (gmx_strncasecmp("bondedtypes",header
,5)==0) {
346 get_a_line(in
,line
,STRLEN
);
347 if ((nparam
=sscanf(line
,"%d %d %d %d %d %d %d %d",
348 &header_settings
->rb
[ebtsBONDS
].type
,&header_settings
->rb
[ebtsANGLES
].type
,
349 &header_settings
->rb
[ebtsPDIHS
].type
,&header_settings
->rb
[ebtsIDIHS
].type
,
350 &dum1
,&header_settings
->nrexcl
,&dum2
,&dum3
)) < 4 )
352 gmx_fatal(FARGS
,"need 4 to 8 parameters in the header of .rtp file %s at line:\n%s\n", rrdb
, line
);
354 header_settings
->bKeepAllGeneratedDihedrals
= (dum1
!= 0);
355 header_settings
->bGenerateHH14Interactions
= (dum2
!= 0);
356 header_settings
->bRemoveDihedralIfWithImproper
= (dum3
!= 0);
357 get_a_line(in
,line
,STRLEN
);
359 fprintf(stderr
,"Using default: not generating all possible dihedrals\n");
360 header_settings
->bKeepAllGeneratedDihedrals
= FALSE
;
363 fprintf(stderr
,"Using default: excluding 3 bonded neighbors\n");
364 header_settings
->nrexcl
= 3;
367 fprintf(stderr
,"Using default: generating 1,4 H--H interactions\n");
368 header_settings
->bGenerateHH14Interactions
= TRUE
;
371 fprintf(stderr
,"Using default: removing proper dihedrals found on the same bond as a proper dihedral\n");
372 header_settings
->bRemoveDihedralIfWithImproper
= TRUE
;
376 "Reading .rtp file without '[ bondedtypes ]' directive,\n"
377 "Will proceed as if the entry was:\n");
378 print_resall_header(stderr
, header_settings
);
380 /* We don't know the current size of rrtp, but simply realloc immediately */
385 if (nrtp
>= maxrtp
) {
390 /* Initialise rtp entry structure */
391 rrtp
[nrtp
] = *header_settings
;
392 if (!get_header(line
,header
))
393 gmx_fatal(FARGS
,"in .rtp file at line:\n%s\n",line
);
394 rrtp
[nrtp
].resname
= strdup(header
);
395 rrtp
[nrtp
].filebase
= strdup(filebase
);
397 get_a_line(in
,line
,STRLEN
);
401 if (!get_header(line
,header
)) {
406 /* header is an bonded directive */
407 bError
= !read_bondeds(bt
,in
,line
,&rrtp
[nrtp
]);
408 } else if (gmx_strncasecmp("atoms",header
,5) == 0) {
409 /* header is the atoms directive */
410 bError
= !read_atoms(in
,line
,&(rrtp
[nrtp
]),tab
,atype
);
412 /* else header must be a residue name */
417 gmx_fatal(FARGS
,"in .rtp file in residue %s at line:\n%s\n",
418 rrtp
[nrtp
].resname
,line
);
419 } while (!feof(in
) && !bNextResidue
);
421 if (rrtp
[nrtp
].natom
== 0)
422 gmx_fatal(FARGS
,"No atoms found in .rtp file in residue %s\n",
425 fprintf(debug
,"%3d %5s %5d",
426 nrtp
+1,rrtp
[nrtp
].resname
,rrtp
[nrtp
].natom
);
427 for(i
=0; i
<ebtsNR
; i
++)
428 fprintf(debug
," %10d",rrtp
[nrtp
].rb
[i
].nb
);
433 for(i
=0; i
<nrtp
; i
++) {
434 if (gmx_strcasecmp(rrtp
[i
].resname
,rrtp
[nrtp
].resname
) == 0) {
439 if (firstrtp
== -1) {
441 fprintf(stderr
,"\rResidue %d",nrtp
);
443 if (firstrtp
>= *nrtpptr
)
445 gmx_fatal(FARGS
,"Found a second entry for '%s' in '%s'",
446 rrtp
[nrtp
].resname
,rrdb
);
448 if (bAllowOverrideRTP
)
450 fprintf(stderr
,"\n");
451 fprintf(stderr
,"Found another rtp entry for '%s' in '%s', ignoring this entry and keeping the one from '%s.rtp'\n",
452 rrtp
[nrtp
].resname
,rrdb
,rrtp
[firstrtp
].filebase
);
453 /* We should free all the data for this entry.
454 * The current code gives a lot of dangling pointers.
456 clear_t_restp(&rrtp
[nrtp
]);
460 gmx_fatal(FARGS
,"Found rtp entries for '%s' in both '%s' and '%s'. If you want the first definition to override the second one, set the -rtpo option of pdb2gmx.",rrtp
[nrtp
].resname
,rrtp
[firstrtp
].filebase
,rrdb
);
465 /* give back unused memory */
468 fprintf(stderr
,"\nSorting it all out...\n");
469 qsort(rrtp
,nrtp
,(size_t)sizeof(rrtp
[0]),comprtp
);
471 check_rtp(nrtp
,rrtp
,rrdb
);
477 /************************************************************
481 ***********************************************************/
482 static gmx_bool
is_sign(char c
)
484 return (c
== '+' || c
== '-');
487 /* Compares if the strins match except for a sign at the end
488 * of (only) one of the two.
490 static int neq_str_sign(const char *a1
,const char *a2
)
494 l1
= (int)strlen(a1
);
495 l2
= (int)strlen(a2
);
499 ((l1
== l2
+1 && is_sign(a1
[l1
-1])) ||
500 (l2
== l1
+1 && is_sign(a2
[l2
-1]))) &&
501 gmx_strncasecmp(a1
,a2
,lm
) == 0)
511 char *search_rtp(const char *key
,int nrtp
,t_restp rtp
[])
513 int i
,n
,nbest
,best
,besti
;
514 char bestbuf
[STRLEN
];
518 /* We want to match at least one character */
520 for(i
=0; (i
<nrtp
); i
++)
522 if (gmx_strcasecmp(key
,rtp
[i
].resname
) == 0)
530 /* Allow a mismatch of at most a sign character (with warning) */
531 n
= neq_str_sign(key
,rtp
[i
].resname
);
533 n
+1 >= (int)strlen(key
) &&
534 n
+1 >= (int)strlen(rtp
[i
].resname
))
540 strcpy(bestbuf
,rtp
[besti
].resname
);
544 strcat(bestbuf
," or ");
545 strcat(bestbuf
,rtp
[i
].resname
);
560 gmx_fatal(FARGS
,"Residue '%s' not found in residue topology database, looks a bit like %s",key
,bestbuf
);
562 else if (besti
== -1)
564 gmx_fatal(FARGS
,"Residue '%s' not found in residue topology database",key
);
566 if (gmx_strcasecmp(rtp
[besti
].resname
,key
) != 0)
569 "\nWARNING: '%s' not found in residue topology database, "
570 "trying to use '%s'\n\n", key
, rtp
[besti
].resname
);
573 return rtp
[besti
].resname
;
576 t_restp
*get_restp(const char *rtpname
,int nrtp
,t_restp rtp
[])
581 while (i
< nrtp
&& gmx_strcasecmp(rtpname
,rtp
[i
].resname
) != 0)
587 /* This should never happen, since search_rtp should have been called
588 * before calling get_restp.
590 gmx_fatal(FARGS
,"Residue type '%s' not found in residue topology database",rtpname
);