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,2017,2018, 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.
46 #include "gromacs/math/utilities.h"
47 #include "gromacs/topology/residuetypes.h"
48 #include "gromacs/utility/cstringutil.h"
49 #include "gromacs/utility/fatalerror.h"
50 #include "gromacs/utility/futil.h"
51 #include "gromacs/utility/pleasecite.h"
52 #include "gromacs/utility/programcontext.h"
53 #include "gromacs/utility/smalloc.h"
54 #include "gromacs/utility/strdb.h"
67 typedef struct gmx_atomprop {
68 gmx_bool bWarned, bWarnVDW;
69 aprop_t prop[epropNR];
70 gmx_residuetype_t *restype;
75 /* NOTFOUND should be smallest, others larger in increasing priority */
77 NOTFOUND = -4, WILDCARD, WILDPROT, PROTEIN
80 /* return number of matching characters,
81 or NOTFOUND if not at least all characters in char *database match */
82 static int dbcmp_len(const char *search, const char *database)
87 while (search[i] && database[i] && (search[i] == database[i]) )
99 static int get_prop_index(aprop_t *ap, gmx_residuetype_t *restype,
100 char *resnm, char *atomnm,
105 long int malen, mrlen;
106 gmx_bool bProtein, bProtWild;
108 bProtein = gmx_residuetype_is_protein(restype, resnm);
109 bProtWild = (strcmp(resnm, "AAA") == 0);
112 for (i = 0; (i < ap->nprop); i++)
114 rlen = dbcmp_len(resnm, ap->resnm[i]);
115 if (rlen == NOTFOUND)
117 if ( (strcmp(ap->resnm[i], "*") == 0) ||
118 (strcmp(ap->resnm[i], "???") == 0) )
122 else if (strcmp(ap->resnm[i], "AAA") == 0)
127 alen = dbcmp_len(atomnm, ap->atomnm[i]);
128 if ( (alen > NOTFOUND) && (rlen > NOTFOUND))
130 if ( ( (alen > malen) && (rlen >= mrlen)) ||
131 ( (rlen > mrlen) && (alen >= malen) ) )
140 *bExact = ((malen == static_cast<long int>(strlen(atomnm))) &&
141 ((mrlen == static_cast<long int>(strlen(resnm))) ||
142 ((mrlen == WILDPROT) && bProtWild) ||
143 ((mrlen == WILDCARD) && !bProtein && !bProtWild)));
147 fprintf(debug, "searching residue: %4s atom: %4s\n", resnm, atomnm);
150 fprintf(debug, " not successful\n");
154 fprintf(debug, " match: %4s %4s\n", ap->resnm[j], ap->atomnm[j]);
160 static void add_prop(aprop_t *ap, gmx_residuetype_t *restype,
161 char *resnm, char *atomnm,
167 j = get_prop_index(ap, restype, resnm, atomnm, &bExact);
171 if (ap->nprop >= ap->maxprop)
174 srenew(ap->resnm, ap->maxprop);
175 srenew(ap->atomnm, ap->maxprop);
176 srenew(ap->value, ap->maxprop);
177 srenew(ap->bAvail, ap->maxprop);
178 for (i = ap->nprop; (i < ap->maxprop); i++)
180 ap->atomnm[i] = nullptr;
181 ap->resnm[i] = nullptr;
183 ap->bAvail[i] = FALSE;
186 ap->atomnm[ap->nprop] = gmx_strdup(atomnm);
187 ap->resnm[ap->nprop] = gmx_strdup(resnm);
193 if (ap->value[j] == p)
195 fprintf(stderr, "Warning double identical entries for %s %s %g on line %d in file %s\n",
196 resnm, atomnm, p, line, ap->db);
200 fprintf(stderr, "Warning double different entries %s %s %g and %g on line %d in file %s\n"
201 "Using last entry (%g)\n",
202 resnm, atomnm, p, ap->value[j], line, ap->db, p);
208 ap->bAvail[j] = TRUE;
213 static void read_prop(gmx_atomprop_t aps, int eprop, double factor)
215 gmx_atomprop *ap2 = static_cast<gmx_atomprop*>(aps);
217 char line[STRLEN], resnm[32], atomnm[32];
222 ap = &ap2->prop[eprop];
224 fp = libopen(ap->db);
226 while (get_a_line(fp, line, STRLEN))
229 if (sscanf(line, "%31s %31s %20lf", resnm, atomnm, &pp) == 3)
232 add_prop(ap, aps->restype, resnm, atomnm, pp, line_no);
236 fprintf(stderr, "WARNING: Error in file %s at line %d ignored\n",
241 /* for libraries we can use the low-level close routines */
247 static void set_prop(gmx_atomprop_t aps, int eprop)
249 gmx_atomprop *ap2 = static_cast<gmx_atomprop*>(aps);
250 const char *fns[epropNR] = { "atommass.dat", "vdwradii.dat", "dgsolv.dat", "electroneg.dat", "elements.dat" };
251 double fac[epropNR] = { 1.0, 1.0, 418.4, 1.0, 1.0 };
252 double def[epropNR] = { 12.011, 0.14, 0.0, 2.2, -1 };
255 ap = &ap2->prop[eprop];
258 ap->db = gmx_strdup(fns[eprop]);
259 ap->def = def[eprop];
260 read_prop(aps, eprop, fac[eprop]);
264 fprintf(debug, "Entries in %s: %d\n", ap->db, ap->nprop);
267 if ( ( (!aps->bWarned) && (eprop == epropMass) ) || (eprop == epropVDW))
270 "WARNING: Masses and atomic (Van der Waals) radii will be guessed\n"
271 " based on residue and atom names, since they could not be\n"
272 " definitively assigned from the information in your input\n"
273 " files. These guessed numbers might deviate from the mass\n"
274 " and radius of the atom type. Please check the output\n"
275 " files if necessary.\n\n");
281 gmx_atomprop_t gmx_atomprop_init(void)
287 gmx_residuetype_init(&aps->restype);
288 aps->bWarned = FALSE;
289 aps->bWarnVDW = FALSE;
291 return static_cast<gmx_atomprop_t>(aps);
294 static void destroy_prop(aprop_t *ap)
302 for (i = 0; i < ap->nprop; i++)
304 sfree(ap->atomnm[i]);
314 void gmx_atomprop_destroy(gmx_atomprop_t aps)
316 gmx_atomprop *ap = static_cast<gmx_atomprop*>(aps);
321 printf("\nWARNING: gmx_atomprop_destroy called with a NULL pointer\n\n");
325 for (p = 0; p < epropNR; p++)
327 destroy_prop(&ap->prop[p]);
330 gmx_residuetype_destroy(ap->restype);
335 static void vdw_warning(FILE *fp)
339 fprintf(fp, "NOTE: From version 5.0 %s uses the Van der Waals radii\n",
340 gmx::getProgramContext().displayName());
341 fprintf(fp, "from the source below. This means the results may be different\n");
342 fprintf(fp, "compared to previous GROMACS versions.\n");
343 please_cite(fp, "Bondi1964a");
347 gmx_bool gmx_atomprop_query(gmx_atomprop_t aps,
348 int eprop, const char *resnm, const char *atomnm,
351 gmx_atomprop *ap = static_cast<gmx_atomprop*>(aps);
354 char atomname[MAXQ], resname[MAXQ];
357 set_prop(aps, eprop);
358 if ((strlen(atomnm) > MAXQ-1) || (strlen(resnm) > MAXQ-1))
362 fprintf(debug, "WARNING: will only compare first %d characters\n",
366 if (isdigit(atomnm[0]))
369 /* put digit after atomname */
370 for (i = 1; i < MAXQ-1 && atomnm[i] != '\0'; i++)
372 atomname[i-1] = atomnm[i];
374 atomname[i-1] = atomnm[0];
379 strncpy(atomname, atomnm, MAXQ-1);
381 strncpy(resname, resnm, MAXQ-1);
383 j = get_prop_index(&(ap->prop[eprop]), ap->restype, resname,
386 if (eprop == epropVDW && !ap->bWarnVDW)
393 *value = ap->prop[eprop].value[j];
398 *value = ap->prop[eprop].def;
403 char *gmx_atomprop_element(gmx_atomprop_t aps, int atomnumber)
405 gmx_atomprop *ap = static_cast<gmx_atomprop*>(aps);
408 set_prop(aps, epropElement);
409 for (i = 0; (i < ap->prop[epropElement].nprop); i++)
411 if (std::round(ap->prop[epropElement].value[i]) == atomnumber)
413 return ap->prop[epropElement].atomnm[i];
419 int gmx_atomprop_atomnumber(gmx_atomprop_t aps, const char *elem)
421 gmx_atomprop *ap = static_cast<gmx_atomprop*>(aps);
424 set_prop(aps, epropElement);
425 for (i = 0; (i < ap->prop[epropElement].nprop); i++)
427 if (gmx_strcasecmp(ap->prop[epropElement].atomnm[i], elem) == 0)
429 return std::round(ap->prop[epropElement].value[i]);