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