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, 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.
44 #include "gromacs/fileio/gmxfio.h"
45 #include "gromacs/fileio/strdb.h"
46 #include "gromacs/gmxpreprocess/fflibutil.h"
47 #include "gromacs/gmxpreprocess/h_db.h"
48 #include "gromacs/gmxpreprocess/notset.h"
49 #include "gromacs/gmxpreprocess/resall.h"
50 #include "gromacs/gmxpreprocess/toputil.h"
51 #include "gromacs/legacyheaders/typedefs.h"
52 #include "gromacs/utility/cstringutil.h"
53 #include "gromacs/utility/fatalerror.h"
54 #include "gromacs/utility/futil.h"
55 #include "gromacs/utility/smalloc.h"
57 /* use bonded types definitions in hackblock.h */
58 #define ekwRepl ebtsNR+1
59 #define ekwAdd ebtsNR+2
60 #define ekwDel ebtsNR+3
62 const char *kw_names
[ekwNR
] = {
63 "replace", "add", "delete"
66 int find_kw(char *keyw
)
70 for (i
= 0; i
< ebtsNR
; i
++)
72 if (gmx_strcasecmp(btsNames
[i
], keyw
) == 0)
77 for (i
= 0; i
< ekwNR
; i
++)
79 if (gmx_strcasecmp(kw_names
[i
], keyw
) == 0)
81 return ebtsNR
+ 1 + i
;
88 #define FATAL() gmx_fatal(FARGS, "Reading Termini Database: not enough items on line\n%s", line)
90 static void read_atom(char *line
, gmx_bool bAdd
,
91 char **nname
, t_atom
*a
, gpp_atomtype_t atype
, int *cgnr
)
97 /* This code is messy, because of support for different formats:
98 * for replace: [new name] <atom type> <m> <q> [cgnr (old format)]
99 * for add: <atom type> <m> <q> [cgnr]
101 nr
= sscanf(line
, "%s %s %s %s %s", buf
[0], buf
[1], buf
[2], buf
[3], buf
[4]);
103 /* Here there an ambiguity due to the old replace format with cgnr,
104 * which was read for years, but ignored in the rest of the code.
105 * We have to assume that the atom type does not start with a digit
106 * to make a line with 4 entries uniquely interpretable.
108 if (!bAdd
&& nr
== 4 && isdigit(buf
[1][0]))
113 if (nr
< 3 || nr
> 4)
115 gmx_fatal(FARGS
, "Reading Termini Database: expected %d or %d items of atom data in stead of %d on line\n%s", 3, 4, nr
, line
);
122 *nname
= gmx_strdup(buf
[i
++]);
129 a
->type
= get_atomtype_type(buf
[i
++], atype
);
130 sscanf(buf
[i
++], "%lf", &m
);
132 sscanf(buf
[i
++], "%lf", &q
);
136 sscanf(buf
[i
++], "%d", cgnr
);
144 static void print_atom(FILE *out
, t_atom
*a
, gpp_atomtype_t atype
)
146 fprintf(out
, "\t%s\t%g\t%g\n",
147 get_atomtype_name(a
->type
, atype
), a
->m
, a
->q
);
150 static void print_ter_db(const char *ff
, char C
, int nb
, t_hackblock tb
[],
151 gpp_atomtype_t atype
)
154 int i
, j
, k
, bt
, nrepl
, nadd
, ndel
;
157 sprintf(buf
, "%s-%c.tdb", ff
, C
);
158 out
= gmx_fio_fopen(buf
, "w");
160 for (i
= 0; (i
< nb
); i
++)
162 fprintf(out
, "[ %s ]\n", tb
[i
].name
);
168 for (j
= 0; j
< tb
[i
].nhack
; j
++)
170 if (tb
[i
].hack
[j
].oname
!= NULL
&& tb
[i
].hack
[j
].nname
!= NULL
)
174 else if (tb
[i
].hack
[j
].oname
== NULL
&& tb
[i
].hack
[j
].nname
!= NULL
)
178 else if (tb
[i
].hack
[j
].oname
!= NULL
&& tb
[i
].hack
[j
].nname
== NULL
)
182 else if (tb
[i
].hack
[j
].oname
== NULL
&& tb
[i
].hack
[j
].nname
== NULL
)
184 gmx_fatal(FARGS
, "invalid hack (%s) in termini database", tb
[i
].name
);
189 fprintf(out
, "[ %s ]\n", kw_names
[ekwRepl
-ebtsNR
-1]);
190 for (j
= 0; j
< tb
[i
].nhack
; j
++)
192 if (tb
[i
].hack
[j
].oname
!= NULL
&& tb
[i
].hack
[j
].nname
!= NULL
)
194 fprintf(out
, "%s\t", tb
[i
].hack
[j
].oname
);
195 print_atom(out
, tb
[i
].hack
[j
].atom
, atype
);
201 fprintf(out
, "[ %s ]\n", kw_names
[ekwAdd
-ebtsNR
-1]);
202 for (j
= 0; j
< tb
[i
].nhack
; j
++)
204 if (tb
[i
].hack
[j
].oname
== NULL
&& tb
[i
].hack
[j
].nname
!= NULL
)
206 print_ab(out
, &(tb
[i
].hack
[j
]), tb
[i
].hack
[j
].nname
);
207 print_atom(out
, tb
[i
].hack
[j
].atom
, atype
);
213 fprintf(out
, "[ %s ]\n", kw_names
[ekwDel
-ebtsNR
-1]);
214 for (j
= 0; j
< tb
[i
].nhack
; j
++)
216 if (tb
[i
].hack
[j
].oname
!= NULL
&& tb
[i
].hack
[j
].nname
== NULL
)
218 fprintf(out
, "%s\n", tb
[i
].hack
[j
].oname
);
222 for (bt
= 0; bt
< ebtsNR
; bt
++)
226 fprintf(out
, "[ %s ]\n", btsNames
[bt
]);
227 for (j
= 0; j
< tb
[i
].rb
[bt
].nb
; j
++)
229 for (k
= 0; k
< btsNiatoms
[bt
]; k
++)
231 fprintf(out
, "%s%s", k
? "\t" : "", tb
[i
].rb
[bt
].b
[j
].a
[k
]);
233 if (tb
[i
].rb
[bt
].b
[j
].s
)
235 fprintf(out
, "\t%s", tb
[i
].rb
[bt
].b
[j
].s
);
246 static void read_ter_db_file(char *fn
,
247 int *ntbptr
, t_hackblock
**tbptr
,
248 gpp_atomtype_t atype
)
250 char filebase
[STRLEN
], *ptr
;
252 char header
[STRLEN
], buf
[STRLEN
], line
[STRLEN
];
254 int i
, j
, n
, ni
, kwnr
, nb
, maxnb
, nh
;
256 fflib_filename_base(fn
, filebase
, STRLEN
);
257 /* Remove the C/N termini extension */
258 ptr
= strrchr(filebase
, '.');
267 fprintf(debug
, "Opened %s\n", fn
);
274 get_a_line(in
, line
, STRLEN
);
277 if (get_header(line
, header
))
279 /* this is a new block, or a new keyword */
280 kwnr
= find_kw(header
);
285 /* here starts a new block */
291 clear_t_hackblock(&tb
[nb
]);
292 tb
[nb
].name
= gmx_strdup(header
);
293 tb
[nb
].filebase
= gmx_strdup(filebase
);
300 gmx_fatal(FARGS
, "reading termini database: "
301 "directive expected before line:\n%s\n"
302 "This might be a file in an old format.", line
);
304 /* this is not a header, so it must be data */
307 /* this is a hack: add/rename/delete atoms */
308 /* make space for hacks */
309 if (tb
[nb
].nhack
>= tb
[nb
].maxhack
)
311 tb
[nb
].maxhack
+= 10;
312 srenew(tb
[nb
].hack
, tb
[nb
].maxhack
);
315 clear_t_hack(&(tb
[nb
].hack
[nh
]));
316 for (i
= 0; i
< 4; i
++)
318 tb
[nb
].hack
[nh
].a
[i
] = NULL
;
324 if (kwnr
== ekwRepl
|| kwnr
== ekwDel
)
326 if (sscanf(line
, "%s%n", buf
, &n
) != 1)
328 gmx_fatal(FARGS
, "Reading Termini Database '%s': "
329 "expected atom name on line\n%s", fn
, line
);
331 tb
[nb
].hack
[nh
].oname
= gmx_strdup(buf
);
332 /* we only replace or delete one atom at a time */
333 tb
[nb
].hack
[nh
].nr
= 1;
335 else if (kwnr
== ekwAdd
)
337 read_ab(line
, fn
, &(tb
[nb
].hack
[nh
]));
338 get_a_line(in
, line
, STRLEN
);
342 gmx_fatal(FARGS
, "unimplemented keyword number %d (%s:%d)",
343 kwnr
, __FILE__
, __LINE__
);
345 if (kwnr
== ekwRepl
|| kwnr
== ekwAdd
)
347 snew(tb
[nb
].hack
[nh
].atom
, 1);
348 read_atom(line
+n
, kwnr
== ekwAdd
,
349 &tb
[nb
].hack
[nh
].nname
, tb
[nb
].hack
[nh
].atom
, atype
,
350 &tb
[nb
].hack
[nh
].cgnr
);
351 if (tb
[nb
].hack
[nh
].nname
== NULL
)
353 if (tb
[nb
].hack
[nh
].oname
!= NULL
)
355 tb
[nb
].hack
[nh
].nname
= gmx_strdup(tb
[nb
].hack
[nh
].oname
);
359 gmx_fatal(FARGS
, "Reading Termini Database '%s': don't know which name the new atom should have on line\n%s", fn
, line
);
364 else if (kwnr
>= 0 && kwnr
< ebtsNR
)
366 /* this is bonded data: bonds, angles, dihedrals or impropers */
367 srenew(tb
[nb
].rb
[kwnr
].b
, tb
[nb
].rb
[kwnr
].nb
+1);
369 for (j
= 0; j
< btsNiatoms
[kwnr
]; j
++)
371 if (sscanf(line
+n
, "%s%n", buf
, &ni
) == 1)
373 tb
[nb
].rb
[kwnr
].b
[tb
[nb
].rb
[kwnr
].nb
].a
[j
] = gmx_strdup(buf
);
377 gmx_fatal(FARGS
, "Reading Termini Database '%s': expected %d atom names (found %d) on line\n%s", fn
, btsNiatoms
[kwnr
], j
-1, line
);
381 for (; j
< MAXATOMLIST
; j
++)
383 tb
[nb
].rb
[kwnr
].b
[tb
[nb
].rb
[kwnr
].nb
].a
[j
] = NULL
;
386 sscanf(line
+n
, "%s", buf
);
387 tb
[nb
].rb
[kwnr
].b
[tb
[nb
].rb
[kwnr
].nb
].s
= gmx_strdup(buf
);
388 tb
[nb
].rb
[kwnr
].nb
++;
392 gmx_fatal(FARGS
, "Reading Termini Database: Expecting a header at line\n"
396 get_a_line(in
, line
, STRLEN
);
407 int read_ter_db(const char *ffdir
, char ter
,
408 t_hackblock
**tbptr
, gpp_atomtype_t atype
)
415 sprintf(ext
, ".%c.tdb", ter
);
417 /* Search for termini database files.
418 * Do not generate an error when none are found.
420 ntdbf
= fflib_search_file_end(ffdir
, ext
, FALSE
, &tdbf
);
423 for (f
= 0; f
< ntdbf
; f
++)
425 read_ter_db_file(tdbf
[f
], &ntb
, tbptr
, atype
);
432 print_ter_db("new", ter
, ntb
, *tbptr
, atype
);
438 t_hackblock
**filter_ter(int nrtp
, t_restp rtp
[],
439 int nb
, t_hackblock tb
[],
444 // TODO Four years later, no force fields have ever used this, so decide status of this feature
445 /* Since some force fields (e.g. OPLS) needs different
446 * atomtypes for different residues there could be a lot
447 * of entries in the databases for specific residues
448 * (e.g. GLY-NH3+, SER-NH3+, ALA-NH3+).
450 * To reduce the database size, we assume that a terminus specifier liker
452 * [ GLY|SER|ALA-NH3+ ]
454 * would cover all of the three residue types above.
455 * Wildcards (*,?) are not OK. Don't worry about e.g. GLU vs. GLUH since
456 * pdb2gmx only uses the first 3 letters when calling this routine.
458 * To automate this, this routines scans a list of termini
459 * for the residue name "resname" and returns an allocated list of
460 * pointers to the termini that could be applied to the
461 * residue in question. The variable pointed to by nret will
462 * contain the number of valid pointers in the list.
463 * Remember to free the list when you are done with it...
467 int i
, j
, n
, none_idx
;
469 char *rtpname_match
, *s
;
472 rtpname_match
= search_rtp(rtpname
, nrtp
, rtp
);
473 restp
= get_restp(rtpname_match
, nrtp
, rtp
);
478 for (i
= 0; i
< nb
; i
++)
484 /* The residue name should appear in a tdb file with the same base name
485 * as the file containing the rtp entry.
486 * This makes termini selection for different molecule types
489 if (gmx_strcasecmp(restp
->filebase
, tb
[i
].filebase
) == 0 &&
490 gmx_strncasecmp(resname
, s
, 3) == 0)
499 /* advance to next |-separated field */
507 while (!found
&& s
!= NULL
);
510 /* All residue-specific termini have been added. We might have to fall
511 * back on generic termini, which are characterized by not having
512 * '-' in the name prior to the last position (which indicates charge).
513 * The [ None ] alternative is special since we don't want that
514 * to be the default, so we put it last in the list we return.
517 for (i
= 0; i
< nb
; i
++)
520 /* The residue name should appear in a tdb file with the same base name
521 * as the file containing the rtp entry.
522 * This makes termini selection for different molecule types
525 if (gmx_strcasecmp(restp
->filebase
, tb
[i
].filebase
) == 0)
527 if (!gmx_strcasecmp("None", s
))
533 /* Time to see if there's a generic terminus that matches.
534 Is there a hyphen? */
535 char *c
= strchr(s
, '-');
537 /* A conjunction hyphen normally indicates a residue-specific
538 terminus, which is named like "GLY-COOH". A generic terminus
539 won't have a hyphen. */
540 bool bFoundAnyHyphen
= (c
!= NULL
);
541 /* '-' as the last character indicates charge, so if that's
542 the only one found e.g. "COO-", then it was not a conjunction
543 hyphen, so this is a generic terminus */
544 bool bOnlyFoundChargeHyphen
= (bFoundAnyHyphen
&&
546 /* Thus, "GLY-COO-" is not recognized as a generic terminus. */
547 bool bFoundGenericTerminus
= !bFoundAnyHyphen
|| bOnlyFoundChargeHyphen
;
548 if (bFoundGenericTerminus
)
550 /* Check that we haven't already added a residue-specific version
553 for (j
= 0; j
< n
&& strstr((*list
[j
]).name
, s
) == NULL
; j
++)
570 list
[n
] = &(tb
[none_idx
]);
579 t_hackblock
*choose_ter(int nb
, t_hackblock
**tb
, const char *title
)
583 printf("%s\n", title
);
584 for (i
= 0; (i
< nb
); i
++)
586 bool bIsZwitterion
= (0 == gmx_wcmatch("*ZWITTERION*", (*tb
[i
]).name
));
587 printf("%2d: %s%s\n", i
, (*tb
[i
]).name
,
588 bIsZwitterion
? " (only use with zwitterions containing exactly one residue)" : "");
592 ret
= fscanf(stdin
, "%d", &sel
);
594 while ((ret
!= 1) || (sel
< 0) || (sel
>= nb
));