src/gromacs/topology/atomprop.cpp

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