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, 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 static void hacksearch_atom(int* ii,
86 gmx::ArrayRef<const std::vector<MoleculePatch>> patches,
98 for (i = 0; (i < pdba->nr) && (pdba->atom[i].resind != resind); i++) {}
99 for (; (i < pdba->nr) && (pdba->atom[i].resind == resind) && (*ii < 0); i++)
102 for (const auto& patch : patches[i])
104 if (patch.nname == name)
114 static std::vector<MoleculePatchDatabase>
115 getMoleculePatchDatabases(const t_atoms* pdba,
116 gmx::ArrayRef<const MoleculePatchDatabase> globalPatches,
118 gmx::ArrayRef<MoleculePatchDatabase* const> ntdb,
119 gmx::ArrayRef<MoleculePatchDatabase* const> ctdb,
120 gmx::ArrayRef<const int> rN,
121 gmx::ArrayRef<const int> rC)
123 std::vector<MoleculePatchDatabase> modBlock(pdba->nres);
125 /* first the termini */
126 for (int i = 0; i < nterpairs; i++)
128 if (ntdb[i] != nullptr)
130 copyModificationBlocks(*ntdb[i], &modBlock[rN[i]]);
132 if (ctdb[i] != nullptr)
134 mergeAtomAndBondModifications(*ctdb[i], &modBlock[rC[i]]);
137 /* then the whole hdb */
138 for (int rnr = 0; rnr < pdba->nres; rnr++)
140 auto ahptr = search_h_db(globalPatches, *pdba->resinfo[rnr].rtp);
141 if (ahptr != globalPatches.end())
143 if (modBlock[rnr].name.empty())
145 modBlock[rnr].name = ahptr->name;
147 mergeAtomModifications(*ahptr, &modBlock[rnr]);
153 static void expand_hackblocks_one(const MoleculePatchDatabase& newPatch,
154 const std::string localAtomName, //NOLINT(performance-unnecessary-value-param)
155 std::vector<MoleculePatch>* globalPatches,
159 /* we'll recursively add atoms to atoms */
161 for (auto& singlePatch : newPatch.hack)
163 /* first check if we're in the N- or C-terminus, then we should ignore
164 all hacks involving atoms from resp. previous or next residue
165 (i.e. which name begins with '-' (N) or '+' (C) */
166 bool bIgnore = false;
167 if (bN) /* N-terminus: ignore '-' */
169 for (int k = 0; k < 4 && !singlePatch.a[k].empty() && !bIgnore; k++)
171 bIgnore = singlePatch.a[k][0] == '-';
174 if (bC) /* C-terminus: ignore '+' */
176 for (int k = 0; k < 4 && !singlePatch.a[k].empty() && !bIgnore; k++)
178 bIgnore = singlePatch.a[k][0] == '+';
181 /* must be either hdb entry (tp>0) or add from tdb (oname==NULL)
182 and first control aton (AI) matches this atom or
183 delete/replace from tdb (oname!=NULL) and oname matches this atom */
186 && (((singlePatch.tp > 0 || singlePatch.oname.empty()) && singlePatch.a[0] == localAtomName)
187 || (singlePatch.oname == localAtomName)))
189 /* now expand all hacks for this atom */
190 for (int k = 0; k < singlePatch.nr; k++)
192 globalPatches->push_back(singlePatch);
193 MoleculePatch* patch = &globalPatches->back();
194 patch->bXSet = false;
195 /* if we're adding (oname==NULL) and don't have a new name (nname)
196 yet, build it from localAtomName */
197 if (patch->nname.empty())
199 if (patch->oname.empty())
201 patch->nname = localAtomName;
202 patch->nname[0] = 'H';
210 "Hack '%s' %d, replacing nname '%s' with '%s' (old name '%s')\n",
211 localAtomName.c_str(),
213 patch->nname.c_str(),
214 singlePatch.nname.c_str(),
215 patch->oname.empty() ? "" : patch->oname.c_str());
217 patch->nname = singlePatch.nname;
220 if (singlePatch.tp == 10 && k == 2)
222 /* This is a water virtual site, not a hydrogen */
223 /* Ugly hardcoded name hack */
224 patch->nname.assign("M");
226 else if (singlePatch.tp == 11 && k >= 2)
228 /* This is a water lone pair, not a hydrogen */
229 /* Ugly hardcoded name hack */
230 patch->nname.assign(gmx::formatString("LP%d", 1 + k - 2));
232 else if (singlePatch.nr > 1)
234 /* adding more than one atom, number them */
235 patch->nname.append(gmx::formatString("%d", 1 + k));
239 /* add hacks to atoms we've just added */
240 if (singlePatch.tp > 0 || singlePatch.oname.empty())
242 for (int k = 0; k < singlePatch.nr; k++)
244 expand_hackblocks_one(
246 globalPatches->at(globalPatches->size() - singlePatch.nr + k).nname,
257 static void expand_hackblocks(const t_atoms* pdba,
258 gmx::ArrayRef<const MoleculePatchDatabase> hb,
259 gmx::ArrayRef<std::vector<MoleculePatch>> patches,
261 gmx::ArrayRef<const int> rN,
262 gmx::ArrayRef<const int> rC)
264 for (int i = 0; i < pdba->nr; i++)
267 for (int j = 0; j < nterpairs && !bN; j++)
269 bN = pdba->atom[i].resind == rN[j];
272 for (int j = 0; j < nterpairs && !bC; j++)
274 bC = pdba->atom[i].resind == rC[j];
277 /* add hacks to this atom */
278 expand_hackblocks_one(hb[pdba->atom[i].resind], *pdba->atomname[i], &patches[i], bN, bC);
282 static int check_atoms_present(const t_atoms* pdba,
283 gmx::ArrayRef<std::vector<MoleculePatch>> patches,
284 gmx::ArrayRef<const int> cyclicBondsIndex)
287 for (int i = 0; i < pdba->nr; i++)
289 int rnr = pdba->atom[i].resind;
290 for (auto patch = patches[i].begin(); patch != patches[i].end(); patch++)
292 switch (patch->type())
294 case MoleculePatchType::Add:
297 /* check if the atom is already present */
298 int k = pdbasearch_atom(patch->nname.c_str(), rnr, pdba, "check", TRUE, cyclicBondsIndex);
301 /* We found the added atom. */
302 patch->bAlreadyPresent = true;
306 patch->bAlreadyPresent = false;
307 /* count how many atoms we'll add */
312 case MoleculePatchType::Delete:
318 case MoleculePatchType::Replace: { break;
320 default: { GMX_THROW(gmx::InternalError("Case not handled"));
328 static void calc_all_pos(const t_atoms* pdba,
329 gmx::ArrayRef<const gmx::RVec> x,
330 gmx::ArrayRef<std::vector<MoleculePatch>> patches,
332 gmx::ArrayRef<const int> cyclicBondsIndex)
336 rvec xa[4]; /* control atoms for calc_h_pos */
337 rvec xh[MAXH]; /* hydrogen positions from calc_h_pos */
341 for (int i = 0; i < pdba->nr; i++)
343 int rnr = pdba->atom[i].resind;
344 for (auto patch = patches[i].begin(); patch != patches[i].end(); patch += patch->nr)
346 GMX_RELEASE_ASSERT(patch < patches[i].end(),
347 "The number of patches in the last patch can not exceed the total "
348 "number of patches");
349 /* check if we're adding: */
350 if (patch->type() == MoleculePatchType::Add && patch->tp > 0)
352 bool bFoundAll = true;
353 for (int m = 0; (m < patch->nctl && bFoundAll); m++)
355 int ia = pdbasearch_atom(patch->a[m].c_str(),
358 bCheckMissing ? "atom" : "check",
363 /* not found in original atoms, might still be in
364 * the patch Instructions (patches) */
365 hacksearch_atom(&ii, &jj, patch->a[m].c_str(), patches, rnr, pdba);
368 copy_rvec(patches[ii][jj].newx, xa[m]);
376 "Atom %s not found in residue %s %d"
378 " while adding hydrogens",
380 *pdba->resinfo[rnr].name,
381 pdba->resinfo[rnr].nr,
382 *pdba->resinfo[rnr].rtp);
388 copy_rvec(x[ia], xa[m]);
393 for (int m = 0; (m < MAXH); m++)
395 for (int d = 0; d < DIM; d++)
407 calc_h_pos(patch->tp, xa, xh, &l);
408 for (int m = 0; m < patch->nr; m++)
410 auto next = patch + m;
411 copy_rvec(xh[m], next->newx);
420 static int add_h_low(t_atoms** initialAtoms,
421 t_atoms** modifiedAtoms,
422 std::vector<gmx::RVec>* xptr,
423 gmx::ArrayRef<const MoleculePatchDatabase> globalPatches,
426 gmx::ArrayRef<MoleculePatchDatabase* const> ntdb,
427 gmx::ArrayRef<MoleculePatchDatabase* const> ctdb,
428 gmx::ArrayRef<const int> rN,
429 gmx::ArrayRef<const int> rC,
430 const bool bCheckMissing,
431 gmx::ArrayRef<const int> cyclicBondsIndex)
434 int newi, natoms, nalreadypresent;
435 std::vector<std::vector<MoleculePatch>> patches;
436 std::vector<gmx::RVec> xn;
438 t_atoms* pdba = *initialAtoms;
440 /* set flags for adding hydrogens (according to hdb) */
444 /* We'll have to do all the hard work */
445 /* first get all the hackblocks for each residue: */
446 std::vector<MoleculePatchDatabase> hb =
447 getMoleculePatchDatabases(pdba, globalPatches, nterpairs, ntdb, ctdb, rN, rC);
449 /* expand the hackblocks to atom level */
450 patches.resize(natoms);
451 expand_hackblocks(pdba, hb, patches, nterpairs, rN, rC);
454 /* Now calc the positions */
455 calc_all_pos(pdba, *xptr, patches, bCheckMissing, cyclicBondsIndex);
457 /* we don't have to add atoms that are already present in initialAtoms,
458 so we will remove them from the patches (MoleculePatch) */
459 nadd = check_atoms_present(pdba, patches, cyclicBondsIndex);
461 /* Copy old atoms, making space for new ones */
464 srenew(*modifiedAtoms, 1);
465 init_t_atoms(*modifiedAtoms, natoms + nadd, FALSE);
466 (*modifiedAtoms)->nres = pdba->nres;
467 srenew((*modifiedAtoms)->resinfo, pdba->nres);
468 std::copy(pdba->resinfo, pdba->resinfo + pdba->nres, (*modifiedAtoms)->resinfo);
475 xn.resize(natoms + nadd);
477 for (int i = 0; (i < natoms); i++)
479 /* check if this atom wasn't scheduled for deletion */
480 if (patches[i].empty() || (!patches[i][0].nname.empty()))
482 if (newi >= natoms + nadd)
484 /*gmx_fatal(FARGS,"Not enough space for adding atoms");*/
486 xn.resize(natoms + nadd);
487 srenew((*modifiedAtoms)->atom, natoms + nadd);
488 srenew((*modifiedAtoms)->atomname, natoms + nadd);
490 copy_atom(pdba, i, (*modifiedAtoms), newi, symtab);
491 copy_rvec((*xptr)[i], xn[newi]);
492 /* process the hacks for this atom */
494 for (auto patch = patches[i].begin(); patch != patches[i].end(); patch++)
496 if (patch->type() == MoleculePatchType::Add) /* add */
499 if (newi >= natoms + nadd)
501 /* gmx_fatal(FARGS,"Not enough space for adding atoms");*/
503 xn.resize(natoms + nadd);
504 srenew((*modifiedAtoms)->atom, natoms + nadd);
505 srenew((*modifiedAtoms)->atomname, natoms + nadd);
507 (*modifiedAtoms)->atom[newi].resind = pdba->atom[i].resind;
509 if (!patch->nname.empty()
510 && (patch->oname.empty() || patch->oname == *(*modifiedAtoms)->atomname[newi]))
513 if (patch->type() == MoleculePatchType::Add && patch->bAlreadyPresent)
515 /* This atom is already present, copy it from the input. */
517 copy_atom(pdba, i + nalreadypresent, (*modifiedAtoms), newi, symtab);
518 copy_rvec((*xptr)[i + nalreadypresent], xn[newi]);
525 "Replacing %d '%s' with (old name '%s') %s\n",
527 ((*modifiedAtoms)->atomname[newi] && *(*modifiedAtoms)->atomname[newi])
528 ? *(*modifiedAtoms)->atomname[newi]
530 patch->oname.empty() ? "" : patch->oname.c_str(),
531 patch->nname.c_str());
533 (*modifiedAtoms)->atomname[newi] = put_symtab(symtab, patch->nname.c_str());
536 copy_rvec(patch->newx, xn[newi]);
543 *(*modifiedAtoms)->atomname[newi],
544 (*modifiedAtoms)->atom[newi].m,
545 (*modifiedAtoms)->atom[newi].q);
550 i += nalreadypresent;
553 (*modifiedAtoms)->nr = newi;
556 *initialAtoms = *modifiedAtoms;
563 int add_h(t_atoms** initialAtoms,
564 t_atoms** localAtoms,
565 std::vector<gmx::RVec>* xptr,
566 gmx::ArrayRef<const MoleculePatchDatabase> globalPatches,
569 gmx::ArrayRef<MoleculePatchDatabase* const> ntdb,
570 gmx::ArrayRef<MoleculePatchDatabase* const> ctdb,
571 gmx::ArrayRef<const int> rN,
572 gmx::ArrayRef<const int> rC,
573 const bool bAllowMissing,
574 gmx::ArrayRef<const int> cyclicBondsIndex)
576 int nold, nnew, niter;
578 /* Here we loop to be able to add atoms to added atoms.
579 * We should not check for missing atoms here.
587 initialAtoms, localAtoms, xptr, globalPatches, symtab, nterpairs, ntdb, ctdb, rN, rC, FALSE, cyclicBondsIndex);
592 "More than 100 iterations of add_h. Maybe you are trying to replace an added "
593 "atom (this is not supported)?");
595 } while (nnew > nold);
599 /* Call add_h_low once more, now only for the missing atoms check */
600 add_h_low(initialAtoms, localAtoms, xptr, globalPatches, symtab, nterpairs, ntdb, ctdb, rN, rC, TRUE, cyclicBondsIndex);