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