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