src/gromacs/gmxlib/atomprop.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  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.
  14
  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.
  19  *
  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.
  26  *
  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.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROningen Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <ctype.h>
  40 #include "sysstuff.h"
  41 #include "smalloc.h"
  42 #include "string2.h"
  43 #include "futil.h"
  44 #include "maths.h"
  45 #include "gmx_fatal.h"
  46 #include "atomprop.h"
  47 #include "maths.h"
  48 #include "macros.h"
  49 #include "index.h"
  50 #include "strdb.h"
  51 #include "copyrite.h"
  52
  53 typedef struct {
  54     gmx_bool    bSet;
  55     int         nprop, maxprop;
  56     char       *db;
  57     double      def;
  58     char      **atomnm;
  59     char      **resnm;
  60     gmx_bool   *bAvail;
  61     real       *value;
  62 } aprop_t;
  63
  64 typedef struct gmx_atomprop {
  65     gmx_bool          bWarned;
  66     aprop_t           prop[epropNR];
  67     gmx_residuetype_t restype;
  68 } gmx_atomprop;
  69
  70
  71
  72 /* NOTFOUND should be smallest, others larger in increasing priority */
  73 enum {
  74     NOTFOUND = -4, WILDCARD, WILDPROT, PROTEIN
  75 };
  76
  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)
  80 {
  81     int i;
  82
  83     i = 0;
  84     while (search[i] && database[i] && (search[i] == database[i]) )
  85     {
  86         i++;
  87     }
  88
  89     if (database[i])
  90     {
  91         i = NOTFOUND;
  92     }
  93     return i;
  94 }
  95
  96 static int get_prop_index(aprop_t *ap, gmx_residuetype_t restype,
  97                           char *resnm, char *atomnm,
  98                           gmx_bool *bExact)
  99 {
 100     int      i, j = NOTFOUND;
 101     long int alen, rlen;
 102     long int malen, mrlen;
 103     gmx_bool bProtein, bProtWild;
 104
 105     bProtein  = gmx_residuetype_is_protein(restype, resnm);
 106     bProtWild = (strcmp(resnm, "AAA") == 0);
 107     malen     = NOTFOUND;
 108     mrlen     = NOTFOUND;
 109     for (i = 0; (i < ap->nprop); i++)
 110     {
 111         rlen = dbcmp_len(resnm, ap->resnm[i]);
 112         if (rlen == NOTFOUND)
 113         {
 114             if ( (strcmp(ap->resnm[i], "*") == 0) ||
 115                  (strcmp(ap->resnm[i], "???") == 0) )
 116             {
 117                 rlen = WILDCARD;
 118             }
 119             else if (strcmp(ap->resnm[i], "AAA") == 0)
 120             {
 121                 rlen = WILDPROT;
 122             }
 123         }
 124         alen = dbcmp_len(atomnm, ap->atomnm[i]);
 125         if ( (alen > NOTFOUND) && (rlen > NOTFOUND))
 126         {
 127             if ( ( (alen > malen) && (rlen >= mrlen)) ||
 128                  ( (rlen > mrlen) && (alen >= malen) ) )
 129             {
 130                 malen = alen;
 131                 mrlen = rlen;
 132                 j     = i;
 133             }
 134         }
 135     }
 136
 137     *bExact = ((malen == (long int)strlen(atomnm)) &&
 138                ((mrlen == (long int)strlen(resnm)) ||
 139                 ((mrlen == WILDPROT) && bProtWild) ||
 140                 ((mrlen == WILDCARD) && !bProtein && !bProtWild)));
 141
 142     if (debug)
 143     {
 144         fprintf(debug, "searching residue: %4s atom: %4s\n", resnm, atomnm);
 145         if (j == NOTFOUND)
 146         {
 147             fprintf(debug, " not successful\n");
 148         }
 149         else
 150         {
 151             fprintf(debug, " match: %4s %4s\n", ap->resnm[j], ap->atomnm[j]);
 152         }
 153     }
 154     return j;
 155 }
 156
 157 static void add_prop(aprop_t *ap, gmx_residuetype_t restype,
 158                      char *resnm, char *atomnm,
 159                      real p, int line)
 160 {
 161     int      i, j;
 162     gmx_bool bExact;
 163
 164     j = get_prop_index(ap, restype, resnm, atomnm, &bExact);
 165
 166     if (!bExact)
 167     {
 168         if (ap->nprop >= ap->maxprop)
 169         {
 170             ap->maxprop += 10;
 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++)
 176             {
 177                 ap->atomnm[i] = NULL;
 178                 ap->resnm[i]  = NULL;
 179                 ap->value[i]  = 0;
 180                 ap->bAvail[i] = FALSE;
 181             }
 182         }
 183         ap->atomnm[ap->nprop] = strdup(atomnm);
 184         ap->resnm[ap->nprop]  = strdup(resnm);
 185         j                     = ap->nprop;
 186         ap->nprop++;
 187     }
 188     if (ap->bAvail[j])
 189     {
 190         if (ap->value[j] == p)
 191         {
 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);
 194         }
 195         else
 196         {
 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);
 200             ap->value[j] = p;
 201         }
 202     }
 203     else
 204     {
 205         ap->bAvail[j] = TRUE;
 206         ap->value[j]  = p;
 207     }
 208 }
 209
 210 static void read_prop(gmx_atomprop_t aps, int eprop, double factor)
 211 {
 212     gmx_atomprop *ap2 = (gmx_atomprop*) aps;
 213     FILE         *fp;
 214     char          line[STRLEN], resnm[32], atomnm[32];
 215     double        pp;
 216     int           line_no;
 217     aprop_t      *ap;
 218
 219     ap = &ap2->prop[eprop];
 220
 221     fp      = libopen(ap->db);
 222     line_no = 0;
 223     while (get_a_line(fp, line, STRLEN))
 224     {
 225         line_no++;
 226         if (sscanf(line, "%s %s %lf", resnm, atomnm, &pp) == 3)
 227         {
 228             pp *= factor;
 229             add_prop(ap, aps->restype, resnm, atomnm, pp, line_no);
 230         }
 231         else
 232         {
 233             fprintf(stderr, "WARNING: Error in file %s at line %d ignored\n",
 234                     ap->db, line_no);
 235         }
 236     }
 237
 238     /* for libraries we can use the low-level close routines */
 239     ffclose(fp);
 240
 241     ap->bSet = TRUE;
 242 }
 243
 244 static void set_prop(gmx_atomprop_t aps, int eprop)
 245 {
 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 };
 250     aprop_t      *ap;
 251
 252     ap = &ap2->prop[eprop];
 253     if (!ap->bSet)
 254     {
 255         ap->db  = strdup(fns[eprop]);
 256         ap->def = def[eprop];
 257         read_prop(aps, eprop, fac[eprop]);
 258
 259         if (debug)
 260         {
 261             fprintf(debug, "Entries in %s: %d\n", ap->db, ap->nprop);
 262         }
 263
 264         if ( ( (!aps->bWarned) && (eprop == epropMass) ) || (eprop == epropVDW))
 265         {
 266             printf("\n"
 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");
 273             aps->bWarned = TRUE;
 274         }
 275     }
 276 }
 277
 278 gmx_atomprop_t gmx_atomprop_init(void)
 279 {
 280     gmx_atomprop *aps;
 281     int           p;
 282
 283     snew(aps, 1);
 284
 285     gmx_residuetype_init(&aps->restype);
 286     aps->bWarned = FALSE;
 287
 288     return (gmx_atomprop_t)aps;
 289 }
 290
 291 static void destroy_prop(aprop_t *ap)
 292 {
 293     int i;
 294
 295     if (ap->bSet)
 296     {
 297         sfree(ap->db);
 298
 299         for (i = 0; i < ap->nprop; i++)
 300         {
 301             sfree(ap->atomnm[i]);
 302             sfree(ap->resnm[i]);
 303         }
 304         sfree(ap->atomnm);
 305         sfree(ap->resnm);
 306         sfree(ap->bAvail);
 307         sfree(ap->value);
 308     }
 309 }
 310
 311 void gmx_atomprop_destroy(gmx_atomprop_t aps)
 312 {
 313     gmx_atomprop *ap = (gmx_atomprop*) aps;
 314     int           p;
 315
 316     if (aps == NULL)
 317     {
 318         printf("\nWARNING: gmx_atomprop_destroy called with a NULL pointer\n\n");
 319         return;
 320     }
 321
 322     for (p = 0; p < epropNR; p++)
 323     {
 324         destroy_prop(&ap->prop[p]);
 325     }
 326
 327     gmx_residuetype_destroy(ap->restype);
 328
 329     sfree(ap);
 330 }
 331
 332 gmx_bool gmx_atomprop_query(gmx_atomprop_t aps,
 333                             int eprop, const char *resnm, const char *atomnm,
 334                             real *value)
 335 {
 336     gmx_atomprop *ap = (gmx_atomprop*) aps;
 337     size_t        i;
 338     int           j;
 339 #define MAXQ 32
 340     char          atomname[MAXQ], resname[MAXQ];
 341     gmx_bool      bExact;
 342
 343     set_prop(aps, eprop);
 344     if ((strlen(atomnm) > MAXQ-1) || (strlen(resnm) > MAXQ-1))
 345     {
 346         if (debug)
 347         {
 348             fprintf(debug, "WARNING: will only compare first %d characters\n",
 349                     MAXQ-1);
 350         }
 351     }
 352     if (isdigit(atomnm[0]))
 353     {
 354         /* put digit after atomname */
 355         for (i = 1; (i < min(MAXQ-1, strlen(atomnm))); i++)
 356         {
 357             atomname[i-1] = atomnm[i];
 358         }
 359         atomname[i-1] = atomnm[0];
 360         atomname[i]   = '\0';
 361     }
 362     else
 363     {
 364         strncpy(atomname, atomnm, MAXQ-1);
 365     }
 366     strncpy(resname, resnm, MAXQ-1);
 367
 368     j = get_prop_index(&(ap->prop[eprop]), ap->restype, resname,
 369                        atomname, &bExact);
 370
 371     if (j >= 0)
 372     {
 373         *value = ap->prop[eprop].value[j];
 374         return TRUE;
 375     }
 376     else
 377     {
 378         *value = ap->prop[eprop].def;
 379         return FALSE;
 380     }
 381 }
 382
 383 char *gmx_atomprop_element(gmx_atomprop_t aps, int atomnumber)
 384 {
 385     gmx_atomprop *ap = (gmx_atomprop*) aps;
 386     int           i;
 387
 388     set_prop(aps, epropElement);
 389     for (i = 0; (i < ap->prop[epropElement].nprop); i++)
 390     {
 391         if (gmx_nint(ap->prop[epropElement].value[i]) == atomnumber)
 392         {
 393             return ap->prop[epropElement].atomnm[i];
 394         }
 395     }
 396     return NULL;
 397 }
 398
 399 int gmx_atomprop_atomnumber(gmx_atomprop_t aps, const char *elem)
 400 {
 401     gmx_atomprop *ap = (gmx_atomprop*) aps;
 402     int           i;
 403
 404     set_prop(aps, epropElement);
 405     for (i = 0; (i < ap->prop[epropElement].nprop); i++)
 406     {
 407         if (gmx_strcasecmp(ap->prop[epropElement].atomnm[i], elem) == 0)
 408         {
 409             return gmx_nint(ap->prop[epropElement].value[i]);
 410         }
 411     }
 412     return NOTSET;
 413 }