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 * GROningen Mixture of Alchemy and Childrens' Stories
43 #include "gromacs/fileio/futil.h"
45 #include "gmx_fatal.h"
64 typedef struct gmx_atomprop {
65 gmx_bool bWarned, bWarnVDW;
66 aprop_t prop[epropNR];
67 gmx_residuetype_t restype;
72 /* NOTFOUND should be smallest, others larger in increasing priority */
74 NOTFOUND = -4, WILDCARD, WILDPROT, PROTEIN
77 /* return number of matching characters,
78 or NOTFOUND if not at least all characters in char *database match */
79 static int dbcmp_len(char *search, char *database)
84 while (search[i] && database[i] && (search[i] == database[i]) )
96 static int get_prop_index(aprop_t *ap, gmx_residuetype_t restype,
97 char *resnm, char *atomnm,
102 long int malen, mrlen;
103 gmx_bool bProtein, bProtWild;
105 bProtein = gmx_residuetype_is_protein(restype, resnm);
106 bProtWild = (strcmp(resnm, "AAA") == 0);
109 for (i = 0; (i < ap->nprop); i++)
111 rlen = dbcmp_len(resnm, ap->resnm[i]);
112 if (rlen == NOTFOUND)
114 if ( (strcmp(ap->resnm[i], "*") == 0) ||
115 (strcmp(ap->resnm[i], "???") == 0) )
119 else if (strcmp(ap->resnm[i], "AAA") == 0)
124 alen = dbcmp_len(atomnm, ap->atomnm[i]);
125 if ( (alen > NOTFOUND) && (rlen > NOTFOUND))
127 if ( ( (alen > malen) && (rlen >= mrlen)) ||
128 ( (rlen > mrlen) && (alen >= malen) ) )
137 *bExact = ((malen == (long int)strlen(atomnm)) &&
138 ((mrlen == (long int)strlen(resnm)) ||
139 ((mrlen == WILDPROT) && bProtWild) ||
140 ((mrlen == WILDCARD) && !bProtein && !bProtWild)));
144 fprintf(debug, "searching residue: %4s atom: %4s\n", resnm, atomnm);
147 fprintf(debug, " not successful\n");
151 fprintf(debug, " match: %4s %4s\n", ap->resnm[j], ap->atomnm[j]);
157 static void add_prop(aprop_t *ap, gmx_residuetype_t restype,
158 char *resnm, char *atomnm,
164 j = get_prop_index(ap, restype, resnm, atomnm, &bExact);
168 if (ap->nprop >= ap->maxprop)
171 srenew(ap->resnm, ap->maxprop);
172 srenew(ap->atomnm, ap->maxprop);
173 srenew(ap->value, ap->maxprop);
174 srenew(ap->bAvail, ap->maxprop);
175 for (i = ap->nprop; (i < ap->maxprop); i++)
177 ap->atomnm[i] = NULL;
180 ap->bAvail[i] = FALSE;
183 ap->atomnm[ap->nprop] = strdup(atomnm);
184 ap->resnm[ap->nprop] = strdup(resnm);
190 if (ap->value[j] == p)
192 fprintf(stderr, "Warning double identical entries for %s %s %g on line %d in file %s\n",
193 resnm, atomnm, p, line, ap->db);
197 fprintf(stderr, "Warning double different entries %s %s %g and %g on line %d in file %s\n"
198 "Using last entry (%g)\n",
199 resnm, atomnm, p, ap->value[j], line, ap->db, p);
205 ap->bAvail[j] = TRUE;
210 static void read_prop(gmx_atomprop_t aps, int eprop, double factor)
212 gmx_atomprop *ap2 = (gmx_atomprop*) aps;
214 char line[STRLEN], resnm[32], atomnm[32];
219 ap = &ap2->prop[eprop];
221 fp = libopen(ap->db);
223 while (get_a_line(fp, line, STRLEN))
226 if (sscanf(line, "%s %s %lf", resnm, atomnm, &pp) == 3)
229 add_prop(ap, aps->restype, resnm, atomnm, pp, line_no);
233 fprintf(stderr, "WARNING: Error in file %s at line %d ignored\n",
238 /* for libraries we can use the low-level close routines */
244 static void set_prop(gmx_atomprop_t aps, int eprop)
246 gmx_atomprop *ap2 = (gmx_atomprop*) aps;
247 const char *fns[epropNR] = { "atommass.dat", "vdwradii.dat", "dgsolv.dat", "electroneg.dat", "elements.dat" };
248 double fac[epropNR] = { 1.0, 1.0, 418.4, 1.0, 1.0 };
249 double def[epropNR] = { 12.011, 0.14, 0.0, 2.2, -1 };
252 ap = &ap2->prop[eprop];
255 ap->db = strdup(fns[eprop]);
256 ap->def = def[eprop];
257 read_prop(aps, eprop, fac[eprop]);
261 fprintf(debug, "Entries in %s: %d\n", ap->db, ap->nprop);
264 if ( ( (!aps->bWarned) && (eprop == epropMass) ) || (eprop == epropVDW))
267 "WARNING: Masses and atomic (Van der Waals) radii will be guessed\n"
268 " based on residue and atom names, since they could not be\n"
269 " definitively assigned from the information in your input\n"
270 " files. These guessed numbers might deviate from the mass\n"
271 " and radius of the atom type. Please check the output\n"
272 " files if necessary.\n\n");
278 gmx_atomprop_t gmx_atomprop_init(void)
285 gmx_residuetype_init(&aps->restype);
286 aps->bWarned = FALSE;
287 aps->bWarnVDW = FALSE;
289 return (gmx_atomprop_t)aps;
292 static void destroy_prop(aprop_t *ap)
300 for (i = 0; i < ap->nprop; i++)
302 sfree(ap->atomnm[i]);
312 void gmx_atomprop_destroy(gmx_atomprop_t aps)
314 gmx_atomprop *ap = (gmx_atomprop*) aps;
319 printf("\nWARNING: gmx_atomprop_destroy called with a NULL pointer\n\n");
323 for (p = 0; p < epropNR; p++)
325 destroy_prop(&ap->prop[p]);
328 gmx_residuetype_destroy(ap->restype);
333 static void vdw_warning(FILE *fp)
337 fprintf(fp, "NOTE: From version 5.0 %s uses the Van der Waals radii\n",
339 fprintf(fp, "from the source below. This means the results may be different\n");
340 fprintf(fp, "compared to previous GROMACS versions.\n");
341 please_cite(fp, "Bondi1964a");
345 gmx_bool gmx_atomprop_query(gmx_atomprop_t aps,
346 int eprop, const char *resnm, const char *atomnm,
349 gmx_atomprop *ap = (gmx_atomprop*) aps;
353 char atomname[MAXQ], resname[MAXQ];
356 set_prop(aps, eprop);
357 if ((strlen(atomnm) > MAXQ-1) || (strlen(resnm) > MAXQ-1))
361 fprintf(debug, "WARNING: will only compare first %d characters\n",
365 if (isdigit(atomnm[0]))
367 /* put digit after atomname */
368 for (i = 1; (i < min(MAXQ-1, strlen(atomnm))); i++)
370 atomname[i-1] = atomnm[i];
372 atomname[i-1] = atomnm[0];
377 strncpy(atomname, atomnm, MAXQ-1);
379 strncpy(resname, resnm, MAXQ-1);
381 j = get_prop_index(&(ap->prop[eprop]), ap->restype, resname,
384 if (eprop == epropVDW && !ap->bWarnVDW)
391 *value = ap->prop[eprop].value[j];
396 *value = ap->prop[eprop].def;
401 char *gmx_atomprop_element(gmx_atomprop_t aps, int atomnumber)
403 gmx_atomprop *ap = (gmx_atomprop*) aps;
406 set_prop(aps, epropElement);
407 for (i = 0; (i < ap->prop[epropElement].nprop); i++)
409 if (gmx_nint(ap->prop[epropElement].value[i]) == atomnumber)
411 return ap->prop[epropElement].atomnm[i];
417 int gmx_atomprop_atomnumber(gmx_atomprop_t aps, const char *elem)
419 gmx_atomprop *ap = (gmx_atomprop*) aps;
422 set_prop(aps, epropElement);
423 for (i = 0; (i < ap->prop[epropElement].nprop); i++)
425 if (gmx_strcasecmp(ap->prop[epropElement].atomnm[i], elem) == 0)
427 return gmx_nint(ap->prop[epropElement].value[i]);