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