2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
48 #include "gromacs/fileio/pdbio.h"
49 #include "gromacs/gmxpreprocess/pdb2top.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/utility/cstringutil.h"
52 #include "gromacs/utility/fatalerror.h"
53 #include "gromacs/utility/smalloc.h"
54 #include "gromacs/utility/strdb.h"
58 std::string firstResidue, secondResidue;
59 std::string firstAtomName, secondAtomName;
60 std::string newFirstResidue, newSecondResidue;
61 int firstBondNumber, secondBondNumber;
71 c = toupper(fgetc(stdin));
72 } while ((c != 'Y') && (c != 'N'));
77 std::vector<SpecialBond> generateSpecialBonds()
79 const char* sbfile = "specbond.dat";
81 std::vector<SpecialBond> specialBonds;
82 char r1buf[32], r2buf[32], a1buf[32], a2buf[32], nr1buf[32], nr2buf[32];
87 int nlines = get_lines(sbfile, &lines);
88 for (int i = 0; (i < nlines); i++)
90 if (sscanf(lines[i], "%s%s%d%s%s%d%lf%s%s", r1buf, a1buf, &nb1, r2buf, a2buf, &nb2, &length,
94 fprintf(stderr, "Invalid line '%s' in %s\n", lines[i], sbfile);
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);
117 fprintf(stderr, "%zu out of %d lines of %s converted successfully\n", specialBonds.size(),
123 static bool is_special(gmx::ArrayRef<const SpecialBond> sb, const char* res, const char* atom)
125 for (const auto& bond : sb)
127 if (((strncmp(bond.firstResidue.c_str(), res, 3) == 0)
128 && (gmx::equalCaseInsensitive(bond.firstAtomName, atom)))
129 || ((strncmp(bond.secondResidue.c_str(), res, 3) == 0)
130 && (gmx::equalCaseInsensitive(bond.secondAtomName, atom))))
138 static bool is_bond(gmx::ArrayRef<const SpecialBond> sb, t_atoms* pdba, int a1, int a2, real d, int* index_sb, bool* bSwap)
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;
146 for (const auto& bond : sb)
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))))
155 if ((0.9 * bond.length < d) && (1.1 * bond.length > d))
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))))
166 if ((0.9 * bond.length < d) && (1.1 * bond.length > d))
176 static void rename_1res(t_atoms* pdba, int resind, const char* newres, bool bVerbose)
180 printf("Using rtp entry %s for %s %d\n", newres, *pdba->resinfo[resind].name,
181 pdba->resinfo[resind].nr);
183 /* this used to free *resname, which messes up the symtab! */
184 snew(pdba->resinfo[resind].rtp, 1);
185 *pdba->resinfo[resind].rtp = gmx_strdup(newres);
188 std::vector<DisulfideBond> makeDisulfideBonds(t_atoms* pdba, rvec x[], bool bInteractive, bool bVerbose)
190 std::vector<SpecialBond> specialBonds = generateSpecialBonds();
191 std::vector<DisulfideBond> bonds;
197 if (!specialBonds.empty())
199 std::vector<int> specialBondResIdxs;
200 std::vector<int> specialBondAtomIdxs;
202 for (int i = 0; (i < pdba->nr); i++)
204 /* Check if this atom is special and if it is not a double atom
205 * in the input that still needs to be removed.
208 if (!specialBondAtomIdxs.empty())
210 prevAtom = specialBondAtomIdxs.back();
212 if (is_special(specialBonds, *pdba->resinfo[pdba->atom[i].resind].name, *pdba->atomname[i])
213 && !(!specialBondAtomIdxs.empty() && pdba->atom[prevAtom].resind == pdba->atom[i].resind
214 && gmx_strcasecmp(*pdba->atomname[prevAtom], *pdba->atomname[i]) == 0))
216 specialBondResIdxs.push_back(pdba->atom[i].resind);
217 specialBondAtomIdxs.push_back(i);
220 /* distance matrix d[nspec][nspec] */
221 int nspec = specialBondAtomIdxs.size();
222 std::vector<std::vector<real>> d(nspec);
223 for (int i = 0; (i < nspec); i++)
226 int ai = specialBondAtomIdxs[i];
227 for (int j = 0; (j < nspec); j++)
229 int aj = specialBondAtomIdxs[j];
230 d[i][j] = std::sqrt(distance2(x[ai], x[aj]));
236 fprintf(stderr, "Special Atom Distance matrix:\n");
237 for (int b = 0; (b < nspec); b += MAXCOL)
239 /* print resname/number column headings */
240 fprintf(stderr, "%8s%8s", "", "");
241 int e = std::min(b + MAXCOL, nspec - 1);
242 for (int i = b; (i < e); i++)
244 sprintf(buf, "%s%d", *pdba->resinfo[pdba->atom[specialBondAtomIdxs[i]].resind].name,
245 pdba->resinfo[specialBondResIdxs[i]].nr);
246 fprintf(stderr, "%8s", buf);
248 fprintf(stderr, "\n");
249 /* print atomname/number column headings */
250 fprintf(stderr, "%8s%8s", "", "");
251 e = std::min(b + MAXCOL, nspec - 1);
252 for (int i = b; (i < e); i++)
254 std::string buf = gmx::formatString("%s%d", *pdba->atomname[specialBondAtomIdxs[i]],
255 specialBondAtomIdxs[i] + 1);
256 fprintf(stderr, "%8s", buf.c_str());
258 fprintf(stderr, "\n");
260 e = std::min(b + MAXCOL, nspec);
261 for (int i = b + 1; (i < nspec); i++)
263 std::string buf = gmx::formatString(
264 "%s%d", *pdba->resinfo[pdba->atom[specialBondAtomIdxs[i]].resind].name,
265 pdba->resinfo[specialBondResIdxs[i]].nr);
266 fprintf(stderr, "%8s", buf.c_str());
267 buf = gmx::formatString("%s%d", *pdba->atomname[specialBondAtomIdxs[i]],
268 specialBondAtomIdxs[i] + 1);
269 fprintf(stderr, "%8s", buf.c_str());
270 int e2 = std::min(i, e);
271 for (int j = b; (j < e2); j++)
273 fprintf(stderr, " %7.3f", d[i][j]);
275 fprintf(stderr, "\n");
280 for (int i = 0; (i < nspec); i++)
282 int ai = specialBondAtomIdxs[i];
283 for (int j = i + 1; (j < nspec); j++)
285 int aj = specialBondAtomIdxs[j];
286 /* Ensure creation of at most nspec special bonds to avoid overflowing bonds[] */
287 if (bonds.size() < specialBondAtomIdxs.size()
288 && is_bond(specialBonds, pdba, ai, aj, d[i][j], &index_sb, &bSwap))
290 fprintf(stderr, "%s %s-%d %s-%d and %s-%d %s-%d%s", bInteractive ? "Link" : "Linking",
291 *pdba->resinfo[pdba->atom[ai].resind].name,
292 pdba->resinfo[specialBondResIdxs[i]].nr, *pdba->atomname[ai], ai + 1,
293 *pdba->resinfo[pdba->atom[aj].resind].name,
294 pdba->resinfo[specialBondResIdxs[j]].nr, *pdba->atomname[aj], aj + 1,
295 bInteractive ? " (y/n) ?" : "...\n");
296 bool bDoit = bInteractive ? yesno() : true;
300 DisulfideBond newBond;
301 /* Store the residue numbers in the bonds array */
302 newBond.firstResidue = specialBondResIdxs[i];
303 newBond.secondResidue = specialBondResIdxs[j];
304 newBond.firstAtom = *pdba->atomname[ai];
305 newBond.secondAtom = *pdba->atomname[aj];
306 bonds.push_back(newBond);
307 /* rename residues */
310 rename_1res(pdba, specialBondResIdxs[i],
311 specialBonds[index_sb].newSecondResidue.c_str(), bVerbose);
312 rename_1res(pdba, specialBondResIdxs[j],
313 specialBonds[index_sb].newFirstResidue.c_str(), bVerbose);
317 rename_1res(pdba, specialBondResIdxs[i],
318 specialBonds[index_sb].newFirstResidue.c_str(), bVerbose);
319 rename_1res(pdba, specialBondResIdxs[j],
320 specialBonds[index_sb].newSecondResidue.c_str(), bVerbose);