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,2019, 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 <algorithm>
  47 #include <memory>
  48
  49 #include "gromacs/math/functions.h"
  50 #include "gromacs/math/utilities.h"
  51 #include "gromacs/topology/residuetypes.h"
  52 #include "gromacs/utility/cstringutil.h"
  53 #include "gromacs/utility/fatalerror.h"
  54 #include "gromacs/utility/futil.h"
  55 #include "gromacs/utility/pleasecite.h"
  56 #include "gromacs/utility/programcontext.h"
  57 #include "gromacs/utility/smalloc.h"
  58 #include "gromacs/utility/strdb.h"
  59
  60 /* NOTFOUND should be smallest, others larger in increasing priority */
  61 enum {
  62     NOTFOUND = -4, WILDCARD, WILDPROT, PROTEIN
  63 };
  64
  65 //! Basic entries in AtomProperty.
  66 struct BaseEntry {
  67     //! Default constructor.
  68     BaseEntry(const std::string &aName, const std::string &rName)
  69         : atomName(aName), residueName(rName), isAvailable(false), value(0.0)
  70     {}
  71     //! Name for atom.
  72     std::string atomName;
  73     //! Name for residue.
  74     std::string residueName;
  75     //! Is property available.
  76     bool        isAvailable;
  77     //! Value set for property.
  78     real        value;
  79 };
  80
  81 //! Conglomeration of atom property entries.
  82 struct AtomProperty {
  83     //! Has property been set.
  84     bool                   isSet = false;
  85     //! Database the property is coming from.
  86     std::string            db;
  87     //! Default value for property.
  88     double                 def = 0.0;
  89     //! Basic entries for properties.
  90     std::vector<BaseEntry> entry;
  91 };
  92
  93 //! Implementation detail type for Atomproperties.
  94 class AtomProperties::Impl
  95 {
  96     public:
  97         //! Should user be warned about error.
  98         bool               bWarned = false;
  99         //! Should user be warned about vdW not found.
 100         bool               bWarnVDW = false;
 101         //! The different atom properties.
 102         AtomProperty       prop[epropNR];
 103         //! The residue types.
 104         ResidueType        restype;
 105 };
 106
 107 /*! \brief
 108  * Find number of matching characters in entry.
 109  *
 110  * If not all characters are matching, return NOTFOUND.
 111  * If the length of the database entry is different from the search,
 112  * also return NOTFOUND.
 113  *
 114  * \param[in] search Entry to compare to database.
 115  * \param[in] database Name of the database entry to compare to.
 116  * \returns Number of matching characters or NOTFOUND.
 117  */
 118 static int compareToDatabase(const std::string &search, const std::string &database)
 119 {
 120     if (database.length() > search.length())
 121     {
 122         return NOTFOUND;
 123     }
 124     size_t matches = 0;
 125     for (size_t i = 0; i < database.length(); i++)
 126     {
 127         if (search[i] == database[i])
 128         {
 129             matches++;
 130         }
 131     }
 132     if (matches == database.length())
 133     {
 134         return matches;
 135     }
 136     else
 137     {
 138         return NOTFOUND;
 139     }
 140 }
 141
 142 /*! \brief
 143  * Finds the index for the property being searched.
 144  *
 145  * \param[in] ap Property to search for.
 146  * \param[in] restype Residuetypes in database.
 147  * \param[in] residueName The name of the residue to look for.
 148  * \param[in] atomName The name of the atom to look for.
 149  * \param[in] bExact Do we have the correct match.
 150  * \returns The index for the property.
 151  */
 152 static int findPropertyIndex(AtomProperty *ap, ResidueType *restype,
 153                              const std::string &residueName, const std::string &atomName,
 154                              gmx_bool *bExact)
 155 {
 156     int      j = NOTFOUND;
 157
 158     bool     bProtein  = restype->namedResidueHasType(residueName, "Protein");
 159     bool     bProtWild = residueName == "AAA";
 160     int      malen     = NOTFOUND;
 161     int      mrlen     = NOTFOUND;
 162     for (size_t i = 0; (i < ap->entry.size()); i++)
 163     {
 164         int rlen = compareToDatabase(residueName, ap->entry[i].residueName);
 165         if (rlen == NOTFOUND)
 166         {
 167             if ( (ap->entry[i].residueName == "*") ||
 168                  (ap->entry[i].residueName == "???")  )
 169             {
 170                 rlen = WILDCARD;
 171             }
 172             else if (ap->entry[i].residueName == "AAA")
 173             {
 174                 rlen = WILDPROT;
 175             }
 176         }
 177         int alen = compareToDatabase(atomName, ap->entry[i].atomName);
 178         if ( (alen > NOTFOUND) && (rlen > NOTFOUND))
 179         {
 180             if ( ( (alen > malen) && (rlen >= mrlen)) ||
 181                  ( (rlen > mrlen) && (alen >= malen) ) )
 182             {
 183                 malen = alen;
 184                 mrlen = rlen;
 185                 j     = i;
 186             }
 187         }
 188     }
 189
 190     *bExact = ((malen == static_cast<long int>(atomName.length())) &&
 191                ((mrlen == static_cast<long int>(residueName.length())) ||
 192                 ((mrlen == WILDPROT) && bProtWild) ||
 193                 ((mrlen == WILDCARD) && !bProtein && !bProtWild)));
 194
 195     if (debug)
 196     {
 197         fprintf(debug, "searching residue: %4s atom: %4s\n", residueName.c_str(),
 198                 atomName.c_str());
 199         if (j == NOTFOUND)
 200         {
 201             fprintf(debug, " not successful\n");
 202         }
 203         else
 204         {
 205             fprintf(debug, " match: %4s %4s\n",
 206                     ap->entry[j].residueName.c_str(), ap->entry[j].atomName.c_str());
 207         }
 208     }
 209     return j;
 210 }
 211
 212 /*! \brief
 213  * Add new property to list.
 214  *
 215  * \param[in] ap Atomproperty to add.
 216  * \param[in] restype Residue type database to use.
 217  * \param[in] residueName Name of the residue.
 218  * \param[in] atomName Name of the atom.
 219  * \param[in] propValue Value of property.
 220  * \param[in] line Where to add property.
 221  */
 222 static void addProperty(AtomProperty *ap, ResidueType *restype,
 223                         const std::string &residueName, const std::string &atomName,
 224                         real propValue, int line)
 225 {
 226     bool bExact;
 227     int  j = findPropertyIndex(ap, restype, residueName, atomName, &bExact);
 228
 229     if (!bExact)
 230     {
 231         ap->entry.emplace_back(BaseEntry(atomName, residueName));
 232
 233         j = ap->entry.size() - 1;
 234     }
 235     if (ap->entry[j].isAvailable)
 236     {
 237         if (ap->entry[j].value == propValue)
 238         {
 239             fprintf(stderr, "Warning double identical entries for %s %s %g on line %d in file %s\n",
 240                     residueName.c_str(), atomName.c_str(), propValue, line, ap->db.c_str());
 241         }
 242         else
 243         {
 244             fprintf(stderr, "Warning double different entries %s %s %g and %g on line %d in file %s\n"
 245                     "Using last entry (%g)\n",
 246                     residueName.c_str(), atomName.c_str(),
 247                     propValue, ap->entry[j].value, line, ap->db.c_str(), propValue);
 248             ap->entry[j].value = propValue;
 249         }
 250     }
 251     else
 252     {
 253         ap->entry[j].isAvailable = TRUE;
 254         ap->entry[j].value       = propValue;
 255     }
 256 }
 257
 258 /*! \brief
 259  * Read property value into structure.
 260  *
 261  * \param[in] ap Atomproperty to be read in.
 262  * \param[in] restype Library of residue types.
 263  * \param[in] factor Scaling factor for property.
 264  */
 265 static void readProperty(AtomProperty *ap, ResidueType *restype, double factor)
 266 {
 267     char          line[STRLEN], resnm[32], atomnm[32];
 268
 269     gmx::FilePtr  fp      = gmx::openLibraryFile(ap->db);
 270     int           line_no = 0;
 271     while (get_a_line(fp.get(), line, STRLEN))
 272     {
 273         line_no++;
 274         double pp;
 275         if (sscanf(line, "%31s %31s %20lf", resnm, atomnm, &pp) == 3)
 276         {
 277             pp *= factor;
 278             addProperty(ap, restype, resnm, atomnm, pp, line_no);
 279         }
 280         else
 281         {
 282             fprintf(stderr, "WARNING: Error in file %s at line %d ignored\n",
 283                     ap->db.c_str(), line_no);
 284         }
 285     }
 286     ap->isSet = TRUE;
 287 }
 288
 289 /*! \brief
 290  * Set value for properties.
 291  *
 292  * \param[in] ap Atomproperty to set.
 293  * \param[in] restype Library of residue types.
 294  * \param[in] eprop Which property to set.
 295  * \param[in] haveBeenWarned If we already set a warning before
 296  * \returns True of warning should be printed.
 297  */
 298 static bool setProperties(AtomProperty *ap, ResidueType *restype, int eprop, bool haveBeenWarned)
 299 {
 300     const char       *fns[epropNR]  = { "atommass.dat", "vdwradii.dat", "dgsolv.dat", "electroneg.dat", "elements.dat" };
 301     double            fac[epropNR]  = { 1.0,    1.0,  418.4, 1.0, 1.0 };
 302     double            def[epropNR]  = { 12.011, 0.14, 0.0, 2.2, -1 };
 303
 304     bool              printWarning = false;
 305     if (!ap->isSet)
 306     {
 307         ap->db  = fns[eprop];
 308         ap->def = def[eprop];
 309         readProperty(ap, restype, fac[eprop]);
 310
 311         if (debug)
 312         {
 313             fprintf(debug, "Entries in %s: %zu\n", ap->db.c_str(), ap->entry.size());
 314         }
 315
 316         if ( (!haveBeenWarned && (eprop == epropMass) ) || (eprop == epropVDW))
 317         {
 318             printWarning = true;
 319         }
 320
 321     }
 322     return printWarning;
 323 }
 324
 325 AtomProperties::AtomProperties()
 326     : impl_(new Impl)
 327 {
 328 }
 329
 330 AtomProperties::~AtomProperties()
 331 {
 332 }
 333
 334 AtomProperty *AtomProperties::prop(int eprop)
 335 {
 336     return &impl_->prop[eprop];
 337 }
 338
 339 ResidueType *AtomProperties::restype()
 340 {
 341     return &impl_->restype;
 342 }
 343
 344 //! Print warning that vdW radii and masses are guessed.
 345 static void printWarning()
 346 {
 347     printf("\n"
 348            "WARNING: Masses and atomic (Van der Waals) radii will be guessed\n"
 349            "         based on residue and atom names, since they could not be\n"
 350            "         definitively assigned from the information in your input\n"
 351            "         files. These guessed numbers might deviate from the mass\n"
 352            "         and radius of the atom type. Please check the output\n"
 353            "         files if necessary.\n\n");
 354 }
 355
 356 static void printvdwWarning(FILE *fp)
 357 {
 358     if (nullptr != fp)
 359     {
 360         fprintf(fp, "NOTE: From version 5.0 %s uses the Van der Waals radii\n",
 361                 gmx::getProgramContext().displayName());
 362         fprintf(fp, "from the source below. This means the results may be different\n");
 363         fprintf(fp, "compared to previous GROMACS versions.\n");
 364         please_cite(fp, "Bondi1964a");
 365     }
 366 }
 367
 368 bool AtomProperties::setAtomProperty(int                eprop,
 369                                      const std::string &residueName,
 370                                      const std::string &atomName,
 371                                      real              *value)
 372 {
 373     int           j;
 374     std::string   tmpAtomName, tmpResidueName;
 375     gmx_bool      bExact;
 376
 377     if (setProperties(prop(eprop), restype(), eprop, impl_->bWarned))
 378     {
 379         printWarning();
 380         impl_->bWarned = true;
 381     }
 382     if (isdigit(atomName[0]))
 383     {
 384         /* put digit after atomname */
 385         tmpAtomName.append(atomName.substr(1));
 386         tmpAtomName.append(1, atomName[0]);
 387     }
 388     else
 389     {
 390         tmpAtomName = atomName;
 391     }
 392     j = findPropertyIndex(&(impl_->prop[eprop]), &impl_->restype, residueName,
 393                           tmpAtomName, &bExact);
 394
 395     if (eprop == epropVDW && !impl_->bWarnVDW)
 396     {
 397         printvdwWarning(stdout);
 398         impl_->bWarnVDW = true;
 399     }
 400     if (j >= 0)
 401     {
 402         *value = impl_->prop[eprop].entry[j].value;
 403         return true;
 404     }
 405     else
 406     {
 407         *value = impl_->prop[eprop].def;
 408         return false;
 409     }
 410 }
 411
 412
 413 const std::string AtomProperties::elementFromAtomNumber(int atomNumber)
 414 {
 415     if (setProperties(prop(epropElement), restype(), epropElement, impl_->bWarned))
 416     {
 417         printWarning();
 418         impl_->bWarned = true;
 419     }
 420     for (const auto &e : prop(epropElement)->entry)
 421     {
 422         if (std::round(e.value) == atomNumber)
 423         {
 424             return e.atomName;
 425         }
 426     }
 427     return "";
 428 }
 429
 430 int AtomProperties::atomNumberFromElement(const char *element)
 431 {
 432     if (setProperties(prop(epropElement), restype(), epropElement, impl_->bWarned))
 433     {
 434         printWarning();
 435         impl_->bWarned = true;
 436     }
 437     for (const auto &e : prop(epropElement)->entry)
 438     {
 439         if (gmx_strcasecmp(e.atomName.c_str(), element) == 0)
 440         {
 441             return gmx::roundToInt(e.value);
 442         }
 443     }
 444     return -1;
 445 }