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