src/gromacs/gmxpreprocess/specbond.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,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,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 "specbond.h"
  41
  42 #include <cctype>
  43 #include <cmath>
  44 #include <cstring>
  45
  46 #include <algorithm>
  47
  48 #include "gromacs/fileio/pdbio.h"
  49 #include "gromacs/gmxpreprocess/pdb2top.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/utility/arrayref.h"
  52 #include "gromacs/utility/cstringutil.h"
  53 #include "gromacs/utility/fatalerror.h"
  54 #include "gromacs/utility/smalloc.h"
  55 #include "gromacs/utility/strdb.h"
  56
  57 struct SpecialBond
  58 {
  59     std::string firstResidue, secondResidue;
  60     std::string firstAtomName, secondAtomName;
  61     std::string newFirstResidue, newSecondResidue;
  62     int         firstBondNumber, secondBondNumber;
  63     real        length;
  64 };
  65
  66 static bool yesno()
  67 {
  68     char c;
  69
  70     do
  71     {
  72         c = toupper(fgetc(stdin));
  73     } while ((c != 'Y') && (c != 'N'));
  74
  75     return (c == 'Y');
  76 }
  77
  78 std::vector<SpecialBond> generateSpecialBonds()
  79 {
  80     const char* sbfile = "specbond.dat";
  81
  82     std::vector<SpecialBond> specialBonds;
  83     char                     r1buf[32], r2buf[32], a1buf[32], a2buf[32], nr1buf[32], nr2buf[32];
  84     double                   length;
  85     int                      nb1, nb2;
  86     char**                   lines;
  87
  88     int nlines = get_lines(sbfile, &lines);
  89     for (int i = 0; (i < nlines); i++)
  90     {
  91         if (sscanf(lines[i], "%s%s%d%s%s%d%lf%s%s", r1buf, a1buf, &nb1, r2buf, a2buf, &nb2, &length, nr1buf, nr2buf)
  92             != 9)
  93         {
  94             fprintf(stderr, "Invalid line '%s' in %s\n", lines[i], sbfile);
  95         }
  96         else
  97         {
  98             SpecialBond newBond;
  99
 100             newBond.firstResidue     = r1buf;
 101             newBond.secondResidue    = r2buf;
 102             newBond.newFirstResidue  = nr1buf;
 103             newBond.newSecondResidue = nr2buf;
 104             newBond.firstAtomName    = a1buf;
 105             newBond.secondAtomName   = a2buf;
 106             newBond.firstBondNumber  = nb1;
 107             newBond.secondBondNumber = nb2;
 108             newBond.length           = length;
 109             specialBonds.push_back(newBond);
 110         }
 111         sfree(lines[i]);
 112     }
 113     if (nlines > 0)
 114     {
 115         sfree(lines);
 116     }
 117     fprintf(stderr, "%zu out of %d lines of %s converted successfully\n", specialBonds.size(), nlines, sbfile);
 118
 119     return specialBonds;
 120 }
 121
 122 static bool is_special(gmx::ArrayRef<const SpecialBond> sb, const char* res, const char* atom)
 123 {
 124     for (const auto& bond : sb)
 125     {
 126         if (((strncmp(bond.firstResidue.c_str(), res, 3) == 0)
 127              && (gmx::equalCaseInsensitive(bond.firstAtomName, atom)))
 128             || ((strncmp(bond.secondResidue.c_str(), res, 3) == 0)
 129                 && (gmx::equalCaseInsensitive(bond.secondAtomName, atom))))
 130         {
 131             return TRUE;
 132         }
 133     }
 134     return FALSE;
 135 }
 136
 137 static bool is_bond(gmx::ArrayRef<const SpecialBond> sb, t_atoms* pdba, int a1, int a2, real d, int* index_sb, bool* bSwap)
 138 {
 139     const char* at1  = *pdba->atomname[a1];
 140     const char* at2  = *pdba->atomname[a2];
 141     const char* res1 = *pdba->resinfo[pdba->atom[a1].resind].name;
 142     const char* res2 = *pdba->resinfo[pdba->atom[a2].resind].name;
 143
 144     int i = 0;
 145     for (const auto& bond : sb)
 146     {
 147         *index_sb = i;
 148         if (((strncmp(bond.firstResidue.c_str(), res1, 3) == 0)
 149              && (gmx::equalCaseInsensitive(bond.firstAtomName, at1))
 150              && (strncmp(bond.secondResidue.c_str(), res2, 3) == 0)
 151              && (gmx::equalCaseInsensitive(bond.secondAtomName, at2))))
 152         {
 153             *bSwap = FALSE;
 154             if ((0.9 * bond.length < d) && (1.1 * bond.length > d))
 155             {
 156                 return TRUE;
 157             }
 158         }
 159         if (((strncmp(bond.firstResidue.c_str(), res2, 3) == 0)
 160              && (gmx::equalCaseInsensitive(bond.firstAtomName, at2))
 161              && (strncmp(bond.secondResidue.c_str(), res1, 3) == 0)
 162              && (gmx::equalCaseInsensitive(bond.secondAtomName, at1))))
 163         {
 164             *bSwap = TRUE;
 165             if ((0.9 * bond.length < d) && (1.1 * bond.length > d))
 166             {
 167                 return TRUE;
 168             }
 169         }
 170         i++;
 171     }
 172     return FALSE;
 173 }
 174
 175 static void rename_1res(t_atoms* pdba, int resind, const char* newres, bool bVerbose)
 176 {
 177     if (bVerbose)
 178     {
 179         printf("Using rtp entry %s for %s %d\n",
 180                newres,
 181                *pdba->resinfo[resind].name,
 182                pdba->resinfo[resind].nr);
 183     }
 184     /* this used to free *resname, which messes up the symtab! */
 185     snew(pdba->resinfo[resind].rtp, 1);
 186     *pdba->resinfo[resind].rtp = gmx_strdup(newres);
 187 }
 188
 189 std::vector<DisulfideBond> makeDisulfideBonds(t_atoms* pdba, rvec x[], bool bInteractive, bool bVerbose)
 190 {
 191     std::vector<SpecialBond>   specialBonds = generateSpecialBonds();
 192     std::vector<DisulfideBond> bonds;
 193     bool                       bSwap;
 194     int                        index_sb;
 195     char                       buf[10];
 196
 197
 198     if (!specialBonds.empty())
 199     {
 200         std::vector<int> specialBondResIdxs;
 201         std::vector<int> specialBondAtomIdxs;
 202
 203         for (int i = 0; (i < pdba->nr); i++)
 204         {
 205             /* Check if this atom is special and if it is not a double atom
 206              * in the input that still needs to be removed.
 207              */
 208             int prevAtom = -1;
 209             if (!specialBondAtomIdxs.empty())
 210             {
 211                 prevAtom = specialBondAtomIdxs.back();
 212             }
 213             if (is_special(specialBonds, *pdba->resinfo[pdba->atom[i].resind].name, *pdba->atomname[i])
 214                 && !(!specialBondAtomIdxs.empty() && pdba->atom[prevAtom].resind == pdba->atom[i].resind
 215                      && gmx_strcasecmp(*pdba->atomname[prevAtom], *pdba->atomname[i]) == 0))
 216             {
 217                 specialBondResIdxs.push_back(pdba->atom[i].resind);
 218                 specialBondAtomIdxs.push_back(i);
 219             }
 220         }
 221         /* distance matrix d[nspec][nspec] */
 222         int                            nspec = specialBondAtomIdxs.size();
 223         std::vector<std::vector<real>> d(nspec);
 224         for (int i = 0; (i < nspec); i++)
 225         {
 226             d[i].resize(nspec);
 227             int ai = specialBondAtomIdxs[i];
 228             for (int j = 0; (j < nspec); j++)
 229             {
 230                 int aj  = specialBondAtomIdxs[j];
 231                 d[i][j] = std::sqrt(distance2(x[ai], x[aj]));
 232             }
 233         }
 234         if (nspec > 1)
 235         {
 236 #define MAXCOL 7
 237             fprintf(stderr, "Special Atom Distance matrix:\n");
 238             for (int b = 0; (b < nspec); b += MAXCOL)
 239             {
 240                 /* print resname/number column headings */
 241                 fprintf(stderr, "%8s%8s", "", "");
 242                 int e = std::min(b + MAXCOL, nspec - 1);
 243                 for (int i = b; (i < e); i++)
 244                 {
 245                     sprintf(buf,
 246                             "%s%d",
 247                             *pdba->resinfo[pdba->atom[specialBondAtomIdxs[i]].resind].name,
 248                             pdba->resinfo[specialBondResIdxs[i]].nr);
 249                     fprintf(stderr, "%8s", buf);
 250                 }
 251                 fprintf(stderr, "\n");
 252                 /* print atomname/number column headings */
 253                 fprintf(stderr, "%8s%8s", "", "");
 254                 e = std::min(b + MAXCOL, nspec - 1);
 255                 for (int i = b; (i < e); i++)
 256                 {
 257                     std::string buf = gmx::formatString(
 258                             "%s%d", *pdba->atomname[specialBondAtomIdxs[i]], specialBondAtomIdxs[i] + 1);
 259                     fprintf(stderr, "%8s", buf.c_str());
 260                 }
 261                 fprintf(stderr, "\n");
 262                 /* print matrix */
 263                 e = std::min(b + MAXCOL, nspec);
 264                 for (int i = b + 1; (i < nspec); i++)
 265                 {
 266                     std::string buf = gmx::formatString(
 267                             "%s%d",
 268                             *pdba->resinfo[pdba->atom[specialBondAtomIdxs[i]].resind].name,
 269                             pdba->resinfo[specialBondResIdxs[i]].nr);
 270                     fprintf(stderr, "%8s", buf.c_str());
 271                     buf = gmx::formatString(
 272                             "%s%d", *pdba->atomname[specialBondAtomIdxs[i]], specialBondAtomIdxs[i] + 1);
 273                     fprintf(stderr, "%8s", buf.c_str());
 274                     int e2 = std::min(i, e);
 275                     for (int j = b; (j < e2); j++)
 276                     {
 277                         fprintf(stderr, " %7.3f", d[i][j]);
 278                     }
 279                     fprintf(stderr, "\n");
 280                 }
 281             }
 282         }
 283
 284         for (int i = 0; (i < nspec); i++)
 285         {
 286             int ai = specialBondAtomIdxs[i];
 287             for (int j = i + 1; (j < nspec); j++)
 288             {
 289                 int aj = specialBondAtomIdxs[j];
 290                 /* Ensure creation of at most nspec special bonds to avoid overflowing bonds[] */
 291                 if (bonds.size() < specialBondAtomIdxs.size()
 292                     && is_bond(specialBonds, pdba, ai, aj, d[i][j], &index_sb, &bSwap))
 293                 {
 294                     fprintf(stderr,
 295                             "%s %s-%d %s-%d and %s-%d %s-%d%s",
 296                             bInteractive ? "Link" : "Linking",
 297                             *pdba->resinfo[pdba->atom[ai].resind].name,
 298                             pdba->resinfo[specialBondResIdxs[i]].nr,
 299                             *pdba->atomname[ai],
 300                             ai + 1,
 301                             *pdba->resinfo[pdba->atom[aj].resind].name,
 302                             pdba->resinfo[specialBondResIdxs[j]].nr,
 303                             *pdba->atomname[aj],
 304                             aj + 1,
 305                             bInteractive ? " (y/n) ?" : "...\n");
 306                     bool bDoit = bInteractive ? yesno() : true;
 307
 308                     if (bDoit)
 309                     {
 310                         DisulfideBond newBond;
 311                         /* Store the residue numbers in the bonds array */
 312                         newBond.firstResidue  = specialBondResIdxs[i];
 313                         newBond.secondResidue = specialBondResIdxs[j];
 314                         newBond.firstAtom     = *pdba->atomname[ai];
 315                         newBond.secondAtom    = *pdba->atomname[aj];
 316                         bonds.push_back(newBond);
 317                         /* rename residues */
 318                         if (bSwap)
 319                         {
 320                             rename_1res(pdba,
 321                                         specialBondResIdxs[i],
 322                                         specialBonds[index_sb].newSecondResidue.c_str(),
 323                                         bVerbose);
 324                             rename_1res(pdba,
 325                                         specialBondResIdxs[j],
 326                                         specialBonds[index_sb].newFirstResidue.c_str(),
 327                                         bVerbose);
 328                         }
 329                         else
 330                         {
 331                             rename_1res(pdba,
 332                                         specialBondResIdxs[i],
 333                                         specialBonds[index_sb].newFirstResidue.c_str(),
 334                                         bVerbose);
 335                             rename_1res(pdba,
 336                                         specialBondResIdxs[j],
 337                                         specialBonds[index_sb].newSecondResidue.c_str(),
 338                                         bVerbose);
 339                         }
 340                     }
 341                 }
 342             }
 343         }
 344     }
 345
 346     return bonds;
 347 }