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-2019,2020,2021, 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.
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.
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.
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.
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.
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.
44 #include "gromacs/fileio/confio.h"
45 #include "gromacs/gmxlib/network.h"
46 #include "gromacs/gmxpreprocess/calch.h"
47 #include "gromacs/gmxpreprocess/h_db.h"
48 #include "gromacs/gmxpreprocess/notset.h"
49 #include "gromacs/gmxpreprocess/pgutil.h"
50 #include "gromacs/gmxpreprocess/ter_db.h"
51 #include "gromacs/math/vec.h"
52 #include "gromacs/topology/atoms.h"
53 #include "gromacs/topology/symtab.h"
54 #include "gromacs/utility/cstringutil.h"
55 #include "gromacs/utility/exceptions.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/futil.h"
58 #include "gromacs/utility/smalloc.h"
60 #include "hackblock.h"
63 static void copy_atom(const t_atoms* atoms1, int a1, t_atoms* atoms2, int a2, t_symtab* symtab)
65 atoms2->atom[a2] = atoms1->atom[a1];
66 atoms2->atomname[a2] = put_symtab(symtab, *atoms1->atomname[a1]);
69 static int pdbasearch_atom(const char* name,
72 const char* searchtype,
74 gmx::ArrayRef<const int> cyclicBondsIndex)
78 for (i = 0; (i < pdba->nr) && (pdba->atom[i].resind != resind); i++) {}
80 return search_atom(name, i, pdba, searchtype, bAllowMissing, cyclicBondsIndex);
83 /*! \brief Return the index of the first atom whose residue index
84 * matches and which has a patch with the given name.
86 * \param[out] ii Index of the first atom in the residue that matches,
87 * -1 if no match occurs.
88 * \param[out] jj Index of the patch that matches,
89 * unchanged if no match occurs.
90 * \param[in] name Name of the desired patch to match
91 * \param[in] patches The patch database to search
92 * \param[in] resind The residue index to match
93 * \param[in] pdba The atoms to work with
95 * \todo The short-circuit logic will be simpler if this returned a
96 * std::pair<int, int> as soon as the first double match is found.
98 static void hacksearch_atom(int* ii,
101 gmx::ArrayRef<const std::vector<MoleculePatch>> patches,
113 for (i = 0; (i < pdba->nr) && (pdba->atom[i].resind != resind); i++) {}
114 for (; (i < pdba->nr) && (pdba->atom[i].resind == resind) && (*ii < 0); i++)
117 for (const auto& patch : patches[i])
119 if (patch.nname == name)
133 static std::vector<MoleculePatchDatabase>
134 getMoleculePatchDatabases(const t_atoms* pdba,
135 gmx::ArrayRef<const MoleculePatchDatabase> globalPatches,
137 gmx::ArrayRef<MoleculePatchDatabase* const> ntdb,
138 gmx::ArrayRef<MoleculePatchDatabase* const> ctdb,
139 gmx::ArrayRef<const int> rN,
140 gmx::ArrayRef<const int> rC)
142 std::vector<MoleculePatchDatabase> modBlock(pdba->nres);
144 /* first the termini */
145 for (int i = 0; i < nterpairs; i++)
147 if (ntdb[i] != nullptr)
149 copyModificationBlocks(*ntdb[i], &modBlock[rN[i]]);
151 if (ctdb[i] != nullptr)
153 mergeAtomAndBondModifications(*ctdb[i], &modBlock[rC[i]]);
156 /* then the whole hdb */
157 for (int rnr = 0; rnr < pdba->nres; rnr++)
159 auto ahptr = search_h_db(globalPatches, *pdba->resinfo[rnr].rtp);
160 if (ahptr != globalPatches.end())
162 if (modBlock[rnr].name.empty())
164 modBlock[rnr].name = ahptr->name;
166 mergeAtomModifications(*ahptr, &modBlock[rnr]);
172 static void expand_hackblocks_one(const MoleculePatchDatabase& newPatch,
173 const std::string localAtomName, //NOLINT(performance-unnecessary-value-param)
174 std::vector<MoleculePatch>* globalPatches,
178 /* we'll recursively add atoms to atoms */
180 for (const auto& singlePatch : newPatch.hack)
182 /* first check if we're in the N- or C-terminus, then we should ignore
183 all hacks involving atoms from resp. previous or next residue
184 (i.e. which name begins with '-' (N) or '+' (C) */
185 bool bIgnore = false;
186 if (bN) /* N-terminus: ignore '-' */
188 for (int k = 0; k < 4 && !singlePatch.a[k].empty() && !bIgnore; k++)
190 bIgnore = singlePatch.a[k][0] == '-';
193 if (bC) /* C-terminus: ignore '+' */
195 for (int k = 0; k < 4 && !singlePatch.a[k].empty() && !bIgnore; k++)
197 bIgnore = singlePatch.a[k][0] == '+';
200 /* must be either hdb entry (tp>0) or add from tdb (oname==NULL)
201 and first control aton (AI) matches this atom or
202 delete/replace from tdb (oname!=NULL) and oname matches this atom */
205 && (((singlePatch.tp > 0 || singlePatch.oname.empty()) && singlePatch.a[0] == localAtomName)
206 || (singlePatch.oname == localAtomName)))
208 /* now expand all hacks for this atom */
209 for (int k = 0; k < singlePatch.nr; k++)
211 globalPatches->push_back(singlePatch);
212 MoleculePatch* patch = &globalPatches->back();
213 patch->bXSet = false;
214 /* if we're adding (oname==NULL) and don't have a new name (nname)
215 yet, build it from localAtomName */
216 if (patch->nname.empty())
218 if (patch->oname.empty())
220 patch->nname = localAtomName;
221 patch->nname[0] = 'H';
229 "Hack '%s' %d, replacing nname '%s' with '%s' (old name '%s')\n",
230 localAtomName.c_str(),
232 patch->nname.c_str(),
233 singlePatch.nname.c_str(),
234 patch->oname.empty() ? "" : patch->oname.c_str());
236 patch->nname = singlePatch.nname;
239 if (singlePatch.tp == 10 && k == 2)
241 /* This is a water virtual site, not a hydrogen */
242 /* Ugly hardcoded name hack to replace 'H' with 'M' */
244 !patch->nname.empty() && patch->nname[0] == 'H',
245 "Water virtual site should be named starting with H at this point");
246 patch->nname[0] = 'M';
248 else if (singlePatch.tp == 11 && k >= 2)
250 /* This is a water lone pair, not a hydrogen */
251 /* Ugly hardcoded name hack */
252 patch->nname.assign(gmx::formatString("LP%d", 1 + k - 2));
254 else if (singlePatch.nr > 1)
256 /* adding more than one atom, number them */
257 patch->nname.append(gmx::formatString("%d", 1 + k));
261 /* add hacks to atoms we've just added */
262 if (singlePatch.tp > 0 || singlePatch.oname.empty())
264 for (int k = 0; k < singlePatch.nr; k++)
266 expand_hackblocks_one(
268 globalPatches->at(globalPatches->size() - singlePatch.nr + k).nname,
279 static void expand_hackblocks(const t_atoms* pdba,
280 gmx::ArrayRef<const MoleculePatchDatabase> hb,
281 gmx::ArrayRef<std::vector<MoleculePatch>> patches,
283 gmx::ArrayRef<const int> rN,
284 gmx::ArrayRef<const int> rC)
286 for (int i = 0; i < pdba->nr; i++)
289 for (int j = 0; j < nterpairs && !bN; j++)
291 bN = pdba->atom[i].resind == rN[j];
294 for (int j = 0; j < nterpairs && !bC; j++)
296 bC = pdba->atom[i].resind == rC[j];
299 /* add hacks to this atom */
300 expand_hackblocks_one(hb[pdba->atom[i].resind], *pdba->atomname[i], &patches[i], bN, bC);
304 static int check_atoms_present(const t_atoms* pdba,
305 gmx::ArrayRef<std::vector<MoleculePatch>> patches,
306 gmx::ArrayRef<const int> cyclicBondsIndex)
309 for (int i = 0; i < pdba->nr; i++)
311 int rnr = pdba->atom[i].resind;
312 for (auto patch = patches[i].begin(); patch != patches[i].end(); patch++)
314 switch (patch->type())
316 case MoleculePatchType::Add:
319 /* check if the atom is already present */
320 int k = pdbasearch_atom(patch->nname.c_str(), rnr, pdba, "check", TRUE, cyclicBondsIndex);
323 /* We found the added atom. */
324 patch->bAlreadyPresent = true;
328 patch->bAlreadyPresent = false;
329 /* count how many atoms we'll add */
334 case MoleculePatchType::Delete:
340 case MoleculePatchType::Replace:
346 GMX_THROW(gmx::InternalError("Case not handled"));
354 static void calc_all_pos(const t_atoms* pdba,
355 gmx::ArrayRef<const gmx::RVec> x,
356 gmx::ArrayRef<std::vector<MoleculePatch>> patches,
358 gmx::ArrayRef<const int> cyclicBondsIndex)
362 rvec xa[4]; /* control atoms for calc_h_pos */
363 rvec xh[MAXH]; /* hydrogen positions from calc_h_pos */
367 for (int i = 0; i < pdba->nr; i++)
369 int rnr = pdba->atom[i].resind;
370 for (auto patch = patches[i].begin(); patch != patches[i].end(); patch += patch->nr)
372 GMX_RELEASE_ASSERT(patch < patches[i].end(),
373 "The number of patches in the last patch can not exceed the total "
374 "number of patches");
375 /* check if we're adding: */
376 if (patch->type() == MoleculePatchType::Add && patch->tp > 0)
378 bool bFoundAll = true;
379 for (int m = 0; (m < patch->nctl && bFoundAll); m++)
381 int ia = pdbasearch_atom(patch->a[m].c_str(),
384 bCheckMissing ? "atom" : "check",
389 /* not found in original atoms, might still be in
390 * the patch Instructions (patches) */
391 hacksearch_atom(&ii, &jj, patch->a[m].c_str(), patches, rnr, pdba);
394 copy_rvec(patches[ii][jj].newx, xa[m]);
402 "Atom %s not found in residue %s %d"
404 " while adding hydrogens",
406 *pdba->resinfo[rnr].name,
407 pdba->resinfo[rnr].nr,
408 *pdba->resinfo[rnr].rtp);
414 copy_rvec(x[ia], xa[m]);
419 for (int m = 0; (m < MAXH); m++)
421 for (int d = 0; d < DIM; d++)
433 calc_h_pos(patch->tp, xa, xh, &l);
434 for (int m = 0; m < patch->nr; m++)
436 auto next = patch + m;
437 copy_rvec(xh[m], next->newx);
446 static int add_h_low(t_atoms** initialAtoms,
447 t_atoms** modifiedAtoms,
448 std::vector<gmx::RVec>* xptr,
449 gmx::ArrayRef<const MoleculePatchDatabase> globalPatches,
452 gmx::ArrayRef<MoleculePatchDatabase* const> ntdb,
453 gmx::ArrayRef<MoleculePatchDatabase* const> ctdb,
454 gmx::ArrayRef<const int> rN,
455 gmx::ArrayRef<const int> rC,
456 const bool bCheckMissing,
457 gmx::ArrayRef<const int> cyclicBondsIndex)
460 int newi, natoms, nalreadypresent;
461 std::vector<std::vector<MoleculePatch>> patches;
462 std::vector<gmx::RVec> xn;
464 t_atoms* pdba = *initialAtoms;
466 /* set flags for adding hydrogens (according to hdb) */
470 /* We'll have to do all the hard work */
471 /* first get all the hackblocks for each residue: */
472 std::vector<MoleculePatchDatabase> hb =
473 getMoleculePatchDatabases(pdba, globalPatches, nterpairs, ntdb, ctdb, rN, rC);
475 /* expand the hackblocks to atom level */
476 patches.resize(natoms);
477 expand_hackblocks(pdba, hb, patches, nterpairs, rN, rC);
480 /* Now calc the positions */
481 calc_all_pos(pdba, *xptr, patches, bCheckMissing, cyclicBondsIndex);
483 /* we don't have to add atoms that are already present in initialAtoms,
484 so we will remove them from the patches (MoleculePatch) */
485 nadd = check_atoms_present(pdba, patches, cyclicBondsIndex);
487 /* Copy old atoms, making space for new ones */
490 srenew(*modifiedAtoms, 1);
491 init_t_atoms(*modifiedAtoms, natoms + nadd, FALSE);
492 (*modifiedAtoms)->nres = pdba->nres;
493 srenew((*modifiedAtoms)->resinfo, pdba->nres);
494 std::copy(pdba->resinfo, pdba->resinfo + pdba->nres, (*modifiedAtoms)->resinfo);
501 xn.resize(natoms + nadd);
503 for (int i = 0; (i < natoms); i++)
505 /* check if this atom wasn't scheduled for deletion */
506 if (patches[i].empty() || (!patches[i][0].nname.empty()))
508 if (newi >= natoms + nadd)
510 /*gmx_fatal(FARGS,"Not enough space for adding atoms");*/
512 xn.resize(natoms + nadd);
513 srenew((*modifiedAtoms)->atom, natoms + nadd);
514 srenew((*modifiedAtoms)->atomname, natoms + nadd);
516 copy_atom(pdba, i, (*modifiedAtoms), newi, symtab);
517 copy_rvec((*xptr)[i], xn[newi]);
518 /* process the hacks for this atom */
520 for (auto patch = patches[i].begin(); patch != patches[i].end(); patch++)
522 if (patch->type() == MoleculePatchType::Add) /* add */
525 if (newi >= natoms + nadd)
527 /* gmx_fatal(FARGS,"Not enough space for adding atoms");*/
529 xn.resize(natoms + nadd);
530 srenew((*modifiedAtoms)->atom, natoms + nadd);
531 srenew((*modifiedAtoms)->atomname, natoms + nadd);
533 (*modifiedAtoms)->atom[newi].resind = pdba->atom[i].resind;
535 if (!patch->nname.empty()
536 && (patch->oname.empty() || patch->oname == *(*modifiedAtoms)->atomname[newi]))
539 if (patch->type() == MoleculePatchType::Add && patch->bAlreadyPresent)
541 /* This atom is already present, copy it from the input. */
543 copy_atom(pdba, i + nalreadypresent, (*modifiedAtoms), newi, symtab);
544 copy_rvec((*xptr)[i + nalreadypresent], xn[newi]);
551 "Replacing %d '%s' with (old name '%s') %s\n",
553 ((*modifiedAtoms)->atomname[newi] && *(*modifiedAtoms)->atomname[newi])
554 ? *(*modifiedAtoms)->atomname[newi]
556 patch->oname.empty() ? "" : patch->oname.c_str(),
557 patch->nname.c_str());
559 (*modifiedAtoms)->atomname[newi] = put_symtab(symtab, patch->nname.c_str());
562 copy_rvec(patch->newx, xn[newi]);
569 *(*modifiedAtoms)->atomname[newi],
570 (*modifiedAtoms)->atom[newi].m,
571 (*modifiedAtoms)->atom[newi].q);
576 i += nalreadypresent;
579 (*modifiedAtoms)->nr = newi;
582 *initialAtoms = *modifiedAtoms;
589 int add_h(t_atoms** initialAtoms,
590 t_atoms** localAtoms,
591 std::vector<gmx::RVec>* xptr,
592 gmx::ArrayRef<const MoleculePatchDatabase> globalPatches,
595 gmx::ArrayRef<MoleculePatchDatabase* const> ntdb,
596 gmx::ArrayRef<MoleculePatchDatabase* const> ctdb,
597 gmx::ArrayRef<const int> rN,
598 gmx::ArrayRef<const int> rC,
599 const bool bAllowMissing,
600 gmx::ArrayRef<const int> cyclicBondsIndex)
602 int nold, nnew, niter;
604 /* Here we loop to be able to add atoms to added atoms.
605 * We should not check for missing atoms here.
613 initialAtoms, localAtoms, xptr, globalPatches, symtab, nterpairs, ntdb, ctdb, rN, rC, FALSE, cyclicBondsIndex);
618 "More than 100 iterations of add_h. Maybe you are trying to replace an added "
619 "atom (this is not supported)?");
621 } while (nnew > nold);
625 /* Call add_h_low once more, now only for the missing atoms check */
626 add_h_low(initialAtoms, localAtoms, xptr, globalPatches, symtab, nterpairs, ntdb, ctdb, rN, rC, TRUE, cyclicBondsIndex);