src/gromacs/gmxlib/atomprop.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013, 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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <ctype.h>
  42 #include "sysstuff.h"
  43 #include "smalloc.h"
  44 #include "string2.h"
  45 #include "gromacs/fileio/futil.h"
  46 #include "maths.h"
  47 #include "gmx_fatal.h"
  48 #include "atomprop.h"
  49 #include "maths.h"
  50 #include "macros.h"
  51 #include "index.h"
  52 #include "strdb.h"
  53 #include "copyrite.h"
  54
  55 typedef struct {
  56     gmx_bool    bSet;
  57     int         nprop, maxprop;
  58     char       *db;
  59     double      def;
  60     char      **atomnm;
  61     char      **resnm;
  62     gmx_bool   *bAvail;
  63     real       *value;
  64 } aprop_t;
  65
  66 typedef struct gmx_atomprop {
  67     gmx_bool          bWarned, bWarnVDW;
  68     aprop_t           prop[epropNR];
  69     gmx_residuetype_t restype;
  70 } gmx_atomprop;
  71
  72
  73
  74 /* NOTFOUND should be smallest, others larger in increasing priority */
  75 enum {
  76     NOTFOUND = -4, WILDCARD, WILDPROT, PROTEIN
  77 };
  78
  79 /* return number of matching characters,
  80    or NOTFOUND if not at least all characters in char *database match */
  81 static int dbcmp_len(char *search, char *database)
  82 {
  83     int i;
  84
  85     i = 0;
  86     while (search[i] && database[i] && (search[i] == database[i]) )
  87     {
  88         i++;
  89     }
  90
  91     if (database[i])
  92     {
  93         i = NOTFOUND;
  94     }
  95     return i;
  96 }
  97
  98 static int get_prop_index(aprop_t *ap, gmx_residuetype_t restype,
  99                           char *resnm, char *atomnm,
 100                           gmx_bool *bExact)
 101 {
 102     int      i, j = NOTFOUND;
 103     long int alen, rlen;
 104     long int malen, mrlen;
 105     gmx_bool bProtein, bProtWild;
 106
 107     bProtein  = gmx_residuetype_is_protein(restype, resnm);
 108     bProtWild = (strcmp(resnm, "AAA") == 0);
 109     malen     = NOTFOUND;
 110     mrlen     = NOTFOUND;
 111     for (i = 0; (i < ap->nprop); i++)
 112     {
 113         rlen = dbcmp_len(resnm, ap->resnm[i]);
 114         if (rlen == NOTFOUND)
 115         {
 116             if ( (strcmp(ap->resnm[i], "*") == 0) ||
 117                  (strcmp(ap->resnm[i], "???") == 0) )
 118             {
 119                 rlen = WILDCARD;
 120             }
 121             else if (strcmp(ap->resnm[i], "AAA") == 0)
 122             {
 123                 rlen = WILDPROT;
 124             }
 125         }
 126         alen = dbcmp_len(atomnm, ap->atomnm[i]);
 127         if ( (alen > NOTFOUND) && (rlen > NOTFOUND))
 128         {
 129             if ( ( (alen > malen) && (rlen >= mrlen)) ||
 130                  ( (rlen > mrlen) && (alen >= malen) ) )
 131             {
 132                 malen = alen;
 133                 mrlen = rlen;
 134                 j     = i;
 135             }
 136         }
 137     }
 138
 139     *bExact = ((malen == (long int)strlen(atomnm)) &&
 140                ((mrlen == (long int)strlen(resnm)) ||
 141                 ((mrlen == WILDPROT) && bProtWild) ||
 142                 ((mrlen == WILDCARD) && !bProtein && !bProtWild)));
 143
 144     if (debug)
 145     {
 146         fprintf(debug, "searching residue: %4s atom: %4s\n", resnm, atomnm);
 147         if (j == NOTFOUND)
 148         {
 149             fprintf(debug, " not successful\n");
 150         }
 151         else
 152         {
 153             fprintf(debug, " match: %4s %4s\n", ap->resnm[j], ap->atomnm[j]);
 154         }
 155     }
 156     return j;
 157 }
 158
 159 static void add_prop(aprop_t *ap, gmx_residuetype_t restype,
 160                      char *resnm, char *atomnm,
 161                      real p, int line)
 162 {
 163     int      i, j;
 164     gmx_bool bExact;
 165
 166     j = get_prop_index(ap, restype, resnm, atomnm, &bExact);
 167
 168     if (!bExact)
 169     {
 170         if (ap->nprop >= ap->maxprop)
 171         {
 172             ap->maxprop += 10;
 173             srenew(ap->resnm, ap->maxprop);
 174             srenew(ap->atomnm, ap->maxprop);
 175             srenew(ap->value, ap->maxprop);
 176             srenew(ap->bAvail, ap->maxprop);
 177             for (i = ap->nprop; (i < ap->maxprop); i++)
 178             {
 179                 ap->atomnm[i] = NULL;
 180                 ap->resnm[i]  = NULL;
 181                 ap->value[i]  = 0;
 182                 ap->bAvail[i] = FALSE;
 183             }
 184         }
 185         ap->atomnm[ap->nprop] = strdup(atomnm);
 186         ap->resnm[ap->nprop]  = strdup(resnm);
 187         j                     = ap->nprop;
 188         ap->nprop++;
 189     }
 190     if (ap->bAvail[j])
 191     {
 192         if (ap->value[j] == p)
 193         {
 194             fprintf(stderr, "Warning double identical entries for %s %s %g on line %d in file %s\n",
 195                     resnm, atomnm, p, line, ap->db);
 196         }
 197         else
 198         {
 199             fprintf(stderr, "Warning double different entries %s %s %g and %g on line %d in file %s\n"
 200                     "Using last entry (%g)\n",
 201                     resnm, atomnm, p, ap->value[j], line, ap->db, p);
 202             ap->value[j] = p;
 203         }
 204     }
 205     else
 206     {
 207         ap->bAvail[j] = TRUE;
 208         ap->value[j]  = p;
 209     }
 210 }
 211
 212 static void read_prop(gmx_atomprop_t aps, int eprop, double factor)
 213 {
 214     gmx_atomprop *ap2 = (gmx_atomprop*) aps;
 215     FILE         *fp;
 216     char          line[STRLEN], resnm[32], atomnm[32];
 217     double        pp;
 218     int           line_no;
 219     aprop_t      *ap;
 220
 221     ap = &ap2->prop[eprop];
 222
 223     fp      = libopen(ap->db);
 224     line_no = 0;
 225     while (get_a_line(fp, line, STRLEN))
 226     {
 227         line_no++;
 228         if (sscanf(line, "%s %s %lf", resnm, atomnm, &pp) == 3)
 229         {
 230             pp *= factor;
 231             add_prop(ap, aps->restype, resnm, atomnm, pp, line_no);
 232         }
 233         else
 234         {
 235             fprintf(stderr, "WARNING: Error in file %s at line %d ignored\n",
 236                     ap->db, line_no);
 237         }
 238     }
 239
 240     /* for libraries we can use the low-level close routines */
 241     ffclose(fp);
 242
 243     ap->bSet = TRUE;
 244 }
 245
 246 static void set_prop(gmx_atomprop_t aps, int eprop)
 247 {
 248     gmx_atomprop *ap2           = (gmx_atomprop*) aps;
 249     const char   *fns[epropNR]  = { "atommass.dat", "vdwradii.dat", "dgsolv.dat", "electroneg.dat", "elements.dat" };
 250     double        fac[epropNR]  = { 1.0,    1.0,  418.4, 1.0, 1.0 };
 251     double        def[epropNR]  = { 12.011, 0.14, 0.0, 2.2, -1 };
 252     aprop_t      *ap;
 253
 254     ap = &ap2->prop[eprop];
 255     if (!ap->bSet)
 256     {
 257         ap->db  = strdup(fns[eprop]);
 258         ap->def = def[eprop];
 259         read_prop(aps, eprop, fac[eprop]);
 260
 261         if (debug)
 262         {
 263             fprintf(debug, "Entries in %s: %d\n", ap->db, ap->nprop);
 264         }
 265
 266         if ( ( (!aps->bWarned) && (eprop == epropMass) ) || (eprop == epropVDW))
 267         {
 268             printf("\n"
 269                    "WARNING: Masses and atomic (Van der Waals) radii will be guessed\n"
 270                    "         based on residue and atom names, since they could not be\n"
 271                    "         definitively assigned from the information in your input\n"
 272                    "         files. These guessed numbers might deviate from the mass\n"
 273                    "         and radius of the atom type. Please check the output\n"
 274                    "         files if necessary.\n\n");
 275             aps->bWarned = TRUE;
 276         }
 277     }
 278 }
 279
 280 gmx_atomprop_t gmx_atomprop_init(void)
 281 {
 282     gmx_atomprop *aps;
 283     int           p;
 284
 285     snew(aps, 1);
 286
 287     gmx_residuetype_init(&aps->restype);
 288     aps->bWarned  = FALSE;
 289     aps->bWarnVDW = FALSE;
 290
 291     return (gmx_atomprop_t)aps;
 292 }
 293
 294 static void destroy_prop(aprop_t *ap)
 295 {
 296     int i;
 297
 298     if (ap->bSet)
 299     {
 300         sfree(ap->db);
 301
 302         for (i = 0; i < ap->nprop; i++)
 303         {
 304             sfree(ap->atomnm[i]);
 305             sfree(ap->resnm[i]);
 306         }
 307         sfree(ap->atomnm);
 308         sfree(ap->resnm);
 309         sfree(ap->bAvail);
 310         sfree(ap->value);
 311     }
 312 }
 313
 314 void gmx_atomprop_destroy(gmx_atomprop_t aps)
 315 {
 316     gmx_atomprop *ap = (gmx_atomprop*) aps;
 317     int           p;
 318
 319     if (aps == NULL)
 320     {
 321         printf("\nWARNING: gmx_atomprop_destroy called with a NULL pointer\n\n");
 322         return;
 323     }
 324
 325     for (p = 0; p < epropNR; p++)
 326     {
 327         destroy_prop(&ap->prop[p]);
 328     }
 329
 330     gmx_residuetype_destroy(ap->restype);
 331
 332     sfree(ap);
 333 }
 334
 335 static void vdw_warning(FILE *fp)
 336 {
 337     if (NULL != fp)
 338     {
 339         fprintf(fp, "NOTE: From version 5.0 %s uses the Van der Waals radii\n",
 340                 ShortProgram());
 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");
 344     }
 345 }
 346
 347 gmx_bool gmx_atomprop_query(gmx_atomprop_t aps,
 348                             int eprop, const char *resnm, const char *atomnm,
 349                             real *value)
 350 {
 351     gmx_atomprop *ap = (gmx_atomprop*) aps;
 352     size_t        i;
 353     int           j;
 354 #define MAXQ 32
 355     char          atomname[MAXQ], resname[MAXQ];
 356     gmx_bool      bExact;
 357
 358     set_prop(aps, eprop);
 359     if ((strlen(atomnm) > MAXQ-1) || (strlen(resnm) > MAXQ-1))
 360     {
 361         if (debug)
 362         {
 363             fprintf(debug, "WARNING: will only compare first %d characters\n",
 364                     MAXQ-1);
 365         }
 366     }
 367     if (isdigit(atomnm[0]))
 368     {
 369         /* put digit after atomname */
 370         for (i = 1; (i < min(MAXQ-1, strlen(atomnm))); i++)
 371         {
 372             atomname[i-1] = atomnm[i];
 373         }
 374         atomname[i-1] = atomnm[0];
 375         atomname[i]   = '\0';
 376     }
 377     else
 378     {
 379         strncpy(atomname, atomnm, MAXQ-1);
 380     }
 381     strncpy(resname, resnm, MAXQ-1);
 382
 383     j = get_prop_index(&(ap->prop[eprop]), ap->restype, resname,
 384                        atomname, &bExact);
 385
 386     if (eprop == epropVDW && !ap->bWarnVDW)
 387     {
 388         vdw_warning(stdout);
 389         ap->bWarnVDW = TRUE;
 390     }
 391     if (j >= 0)
 392     {
 393         *value = ap->prop[eprop].value[j];
 394         return TRUE;
 395     }
 396     else
 397     {
 398         *value = ap->prop[eprop].def;
 399         return FALSE;
 400     }
 401 }
 402
 403 char *gmx_atomprop_element(gmx_atomprop_t aps, int atomnumber)
 404 {
 405     gmx_atomprop *ap = (gmx_atomprop*) aps;
 406     int           i;
 407
 408     set_prop(aps, epropElement);
 409     for (i = 0; (i < ap->prop[epropElement].nprop); i++)
 410     {
 411         if (gmx_nint(ap->prop[epropElement].value[i]) == atomnumber)
 412         {
 413             return ap->prop[epropElement].atomnm[i];
 414         }
 415     }
 416     return NULL;
 417 }
 418
 419 int gmx_atomprop_atomnumber(gmx_atomprop_t aps, const char *elem)
 420 {
 421     gmx_atomprop *ap = (gmx_atomprop*) aps;
 422     int           i;
 423
 424     set_prop(aps, epropElement);
 425     for (i = 0; (i < ap->prop[epropElement].nprop); i++)
 426     {
 427         if (gmx_strcasecmp(ap->prop[epropElement].atomnm[i], elem) == 0)
 428         {
 429             return gmx_nint(ap->prop[epropElement].value[i]);
 430         }
 431     }
 432     return NOTSET;
 433 }