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, 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, int resind, const t_atoms* pdba, const char* searchtype, bool bAllowMissing)
73 for (i = 0; (i < pdba->nr) && (pdba->atom[i].resind != resind); i++) {}
75 return search_atom(name, i, pdba, searchtype, bAllowMissing);
78 static void hacksearch_atom(int* ii,
81 gmx::ArrayRef<const std::vector<MoleculePatch>> patches,
93 for (i = 0; (i < pdba->nr) && (pdba->atom[i].resind != resind); i++) {}
94 for (; (i < pdba->nr) && (pdba->atom[i].resind == resind) && (*ii < 0); i++)
97 for (const auto& patch : patches[i])
99 if (patch.nname == name)
109 static std::vector<MoleculePatchDatabase>
110 getMoleculePatchDatabases(const t_atoms* pdba,
111 gmx::ArrayRef<const MoleculePatchDatabase> globalPatches,
113 gmx::ArrayRef<MoleculePatchDatabase* const> ntdb,
114 gmx::ArrayRef<MoleculePatchDatabase* const> ctdb,
115 gmx::ArrayRef<const int> rN,
116 gmx::ArrayRef<const int> rC)
118 std::vector<MoleculePatchDatabase> modBlock(pdba->nres);
120 /* first the termini */
121 for (int i = 0; i < nterpairs; i++)
123 if (ntdb[i] != nullptr)
125 copyModificationBlocks(*ntdb[i], &modBlock[rN[i]]);
127 if (ctdb[i] != nullptr)
129 mergeAtomAndBondModifications(*ctdb[i], &modBlock[rC[i]]);
132 /* then the whole hdb */
133 for (int rnr = 0; rnr < pdba->nres; rnr++)
135 auto ahptr = search_h_db(globalPatches, *pdba->resinfo[rnr].rtp);
136 if (ahptr != globalPatches.end())
138 if (modBlock[rnr].name.empty())
140 modBlock[rnr].name = ahptr->name;
142 mergeAtomModifications(*ahptr, &modBlock[rnr]);
148 static void expand_hackblocks_one(const MoleculePatchDatabase& newPatch,
149 const std::string localAtomName, //NOLINT(performance-unnecessary-value-param)
150 std::vector<MoleculePatch>* globalPatches,
154 /* we'll recursively add atoms to atoms */
156 for (auto& singlePatch : newPatch.hack)
158 /* first check if we're in the N- or C-terminus, then we should ignore
159 all hacks involving atoms from resp. previous or next residue
160 (i.e. which name begins with '-' (N) or '+' (C) */
161 bool bIgnore = false;
162 if (bN) /* N-terminus: ignore '-' */
164 for (int k = 0; k < 4 && !singlePatch.a[k].empty() && !bIgnore; k++)
166 bIgnore = singlePatch.a[k][0] == '-';
169 if (bC) /* C-terminus: ignore '+' */
171 for (int k = 0; k < 4 && !singlePatch.a[k].empty() && !bIgnore; k++)
173 bIgnore = singlePatch.a[k][0] == '+';
176 /* must be either hdb entry (tp>0) or add from tdb (oname==NULL)
177 and first control aton (AI) matches this atom or
178 delete/replace from tdb (oname!=NULL) and oname matches this atom */
181 && (((singlePatch.tp > 0 || singlePatch.oname.empty()) && singlePatch.a[0] == localAtomName)
182 || (singlePatch.oname == localAtomName)))
184 /* now expand all hacks for this atom */
185 for (int k = 0; k < singlePatch.nr; k++)
187 globalPatches->push_back(singlePatch);
188 MoleculePatch* patch = &globalPatches->back();
189 patch->bXSet = false;
190 /* if we're adding (oname==NULL) and don't have a new name (nname)
191 yet, build it from localAtomName */
192 if (patch->nname.empty())
194 if (patch->oname.empty())
196 patch->nname = localAtomName;
197 patch->nname[0] = 'H';
205 "Hack '%s' %d, replacing nname '%s' with '%s' (old name '%s')\n",
206 localAtomName.c_str(), pos, patch->nname.c_str(),
207 singlePatch.nname.c_str(),
208 patch->oname.empty() ? "" : patch->oname.c_str());
210 patch->nname = singlePatch.nname;
213 if (singlePatch.tp == 10 && k == 2)
215 /* This is a water virtual site, not a hydrogen */
216 /* Ugly hardcoded name hack */
217 patch->nname.assign("M");
219 else if (singlePatch.tp == 11 && k >= 2)
221 /* This is a water lone pair, not a hydrogen */
222 /* Ugly hardcoded name hack */
223 patch->nname.assign(gmx::formatString("LP%d", 1 + k - 2));
225 else if (singlePatch.nr > 1)
227 /* adding more than one atom, number them */
228 patch->nname.append(gmx::formatString("%d", 1 + k));
232 /* add hacks to atoms we've just added */
233 if (singlePatch.tp > 0 || singlePatch.oname.empty())
235 for (int k = 0; k < singlePatch.nr; k++)
237 expand_hackblocks_one(
238 newPatch, globalPatches->at(globalPatches->size() - singlePatch.nr + k).nname,
239 globalPatches, bN, bC);
247 static void expand_hackblocks(const t_atoms* pdba,
248 gmx::ArrayRef<const MoleculePatchDatabase> hb,
249 gmx::ArrayRef<std::vector<MoleculePatch>> patches,
251 gmx::ArrayRef<const int> rN,
252 gmx::ArrayRef<const int> rC)
254 for (int i = 0; i < pdba->nr; i++)
257 for (int j = 0; j < nterpairs && !bN; j++)
259 bN = pdba->atom[i].resind == rN[j];
262 for (int j = 0; j < nterpairs && !bC; j++)
264 bC = pdba->atom[i].resind == rC[j];
267 /* add hacks to this atom */
268 expand_hackblocks_one(hb[pdba->atom[i].resind], *pdba->atomname[i], &patches[i], bN, bC);
272 static int check_atoms_present(const t_atoms* pdba, gmx::ArrayRef<std::vector<MoleculePatch>> patches)
275 for (int i = 0; i < pdba->nr; i++)
277 int rnr = pdba->atom[i].resind;
278 for (auto patch = patches[i].begin(); patch != patches[i].end(); patch++)
280 switch (patch->type())
282 case MoleculePatchType::Add:
285 /* check if the atom is already present */
286 int k = pdbasearch_atom(patch->nname.c_str(), rnr, pdba, "check", TRUE);
289 /* We found the added atom. */
290 patch->bAlreadyPresent = true;
294 patch->bAlreadyPresent = false;
295 /* count how many atoms we'll add */
300 case MoleculePatchType::Delete:
306 case MoleculePatchType::Replace: { break;
308 default: { GMX_THROW(gmx::InternalError("Case not handled"));
316 static void calc_all_pos(const t_atoms* pdba,
317 gmx::ArrayRef<const gmx::RVec> x,
318 gmx::ArrayRef<std::vector<MoleculePatch>> patches,
323 rvec xa[4]; /* control atoms for calc_h_pos */
324 rvec xh[MAXH]; /* hydrogen positions from calc_h_pos */
328 for (int i = 0; i < pdba->nr; i++)
330 int rnr = pdba->atom[i].resind;
331 for (auto patch = patches[i].begin(); patch != patches[i].end(); patch += patch->nr)
333 GMX_RELEASE_ASSERT(patch < patches[i].end(),
334 "The number of patches in the last patch can not exceed the total "
335 "number of patches");
336 /* check if we're adding: */
337 if (patch->type() == MoleculePatchType::Add && patch->tp > 0)
339 bool bFoundAll = true;
340 for (int m = 0; (m < patch->nctl && bFoundAll); m++)
342 int ia = pdbasearch_atom(patch->a[m].c_str(), rnr, pdba,
343 bCheckMissing ? "atom" : "check", !bCheckMissing);
346 /* not found in original atoms, might still be in
347 * the patch Instructions (patches) */
348 hacksearch_atom(&ii, &jj, patch->a[m].c_str(), patches, rnr, pdba);
351 copy_rvec(patches[ii][jj].newx, xa[m]);
359 "Atom %s not found in residue %s %d"
361 " while adding hydrogens",
362 patch->a[m].c_str(), *pdba->resinfo[rnr].name,
363 pdba->resinfo[rnr].nr, *pdba->resinfo[rnr].rtp);
369 copy_rvec(x[ia], xa[m]);
374 for (int m = 0; (m < MAXH); m++)
376 for (int d = 0; d < DIM; d++)
388 calc_h_pos(patch->tp, xa, xh, &l);
389 for (int m = 0; m < patch->nr; m++)
391 auto next = patch + m;
392 copy_rvec(xh[m], next->newx);
401 static int add_h_low(t_atoms** initialAtoms,
402 t_atoms** modifiedAtoms,
403 std::vector<gmx::RVec>* xptr,
404 gmx::ArrayRef<const MoleculePatchDatabase> globalPatches,
407 gmx::ArrayRef<MoleculePatchDatabase* const> ntdb,
408 gmx::ArrayRef<MoleculePatchDatabase* const> ctdb,
409 gmx::ArrayRef<const int> rN,
410 gmx::ArrayRef<const int> rC,
411 const bool bCheckMissing)
414 int newi, natoms, nalreadypresent;
415 std::vector<std::vector<MoleculePatch>> patches;
416 std::vector<gmx::RVec> xn;
418 t_atoms* pdba = *initialAtoms;
420 /* set flags for adding hydrogens (according to hdb) */
424 /* We'll have to do all the hard work */
425 /* first get all the hackblocks for each residue: */
426 std::vector<MoleculePatchDatabase> hb =
427 getMoleculePatchDatabases(pdba, globalPatches, nterpairs, ntdb, ctdb, rN, rC);
429 /* expand the hackblocks to atom level */
430 patches.resize(natoms);
431 expand_hackblocks(pdba, hb, patches, nterpairs, rN, rC);
434 /* Now calc the positions */
435 calc_all_pos(pdba, *xptr, patches, bCheckMissing);
437 /* we don't have to add atoms that are already present in initialAtoms,
438 so we will remove them from the patches (MoleculePatch) */
439 nadd = check_atoms_present(pdba, patches);
441 /* Copy old atoms, making space for new ones */
444 srenew(*modifiedAtoms, 1);
445 init_t_atoms(*modifiedAtoms, natoms + nadd, FALSE);
446 (*modifiedAtoms)->nres = pdba->nres;
447 srenew((*modifiedAtoms)->resinfo, pdba->nres);
448 std::copy(pdba->resinfo, pdba->resinfo + pdba->nres, (*modifiedAtoms)->resinfo);
455 xn.resize(natoms + nadd);
457 for (int i = 0; (i < natoms); i++)
459 /* check if this atom wasn't scheduled for deletion */
460 if (patches[i].empty() || (!patches[i][0].nname.empty()))
462 if (newi >= natoms + nadd)
464 /*gmx_fatal(FARGS,"Not enough space for adding atoms");*/
466 xn.resize(natoms + nadd);
467 srenew((*modifiedAtoms)->atom, natoms + nadd);
468 srenew((*modifiedAtoms)->atomname, natoms + nadd);
470 copy_atom(pdba, i, (*modifiedAtoms), newi, symtab);
471 copy_rvec((*xptr)[i], xn[newi]);
472 /* process the hacks for this atom */
474 for (auto patch = patches[i].begin(); patch != patches[i].end(); patch++)
476 if (patch->type() == MoleculePatchType::Add) /* add */
479 if (newi >= natoms + nadd)
481 /* gmx_fatal(FARGS,"Not enough space for adding atoms");*/
483 xn.resize(natoms + nadd);
484 srenew((*modifiedAtoms)->atom, natoms + nadd);
485 srenew((*modifiedAtoms)->atomname, natoms + nadd);
487 (*modifiedAtoms)->atom[newi].resind = pdba->atom[i].resind;
489 if (!patch->nname.empty()
490 && (patch->oname.empty() || patch->oname == *(*modifiedAtoms)->atomname[newi]))
493 if (patch->type() == MoleculePatchType::Add && patch->bAlreadyPresent)
495 /* This atom is already present, copy it from the input. */
497 copy_atom(pdba, i + nalreadypresent, (*modifiedAtoms), newi, symtab);
498 copy_rvec((*xptr)[i + nalreadypresent], xn[newi]);
504 fprintf(debug, "Replacing %d '%s' with (old name '%s') %s\n", newi,
505 ((*modifiedAtoms)->atomname[newi] && *(*modifiedAtoms)->atomname[newi])
506 ? *(*modifiedAtoms)->atomname[newi]
508 patch->oname.empty() ? "" : patch->oname.c_str(),
509 patch->nname.c_str());
511 (*modifiedAtoms)->atomname[newi] = put_symtab(symtab, patch->nname.c_str());
514 copy_rvec(patch->newx, xn[newi]);
519 fprintf(debug, " %s %g %g", *(*modifiedAtoms)->atomname[newi],
520 (*modifiedAtoms)->atom[newi].m, (*modifiedAtoms)->atom[newi].q);
525 i += nalreadypresent;
528 (*modifiedAtoms)->nr = newi;
531 *initialAtoms = *modifiedAtoms;
538 int add_h(t_atoms** initialAtoms,
539 t_atoms** localAtoms,
540 std::vector<gmx::RVec>* xptr,
541 gmx::ArrayRef<const MoleculePatchDatabase> globalPatches,
544 gmx::ArrayRef<MoleculePatchDatabase* const> ntdb,
545 gmx::ArrayRef<MoleculePatchDatabase* const> ctdb,
546 gmx::ArrayRef<const int> rN,
547 gmx::ArrayRef<const int> rC,
548 const bool bAllowMissing)
550 int nold, nnew, niter;
552 /* Here we loop to be able to add atoms to added atoms.
553 * We should not check for missing atoms here.
560 nnew = add_h_low(initialAtoms, localAtoms, xptr, globalPatches, symtab, nterpairs, ntdb,
561 ctdb, rN, rC, FALSE);
566 "More than 100 iterations of add_h. Maybe you are trying to replace an added "
567 "atom (this is not supported)?");
569 } while (nnew > nold);
573 /* Call add_h_low once more, now only for the missing atoms check */
574 add_h_low(initialAtoms, localAtoms, xptr, globalPatches, symtab, nterpairs, ntdb, ctdb, rN,