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,2017,2018, 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/resall.h"
51 #include "gromacs/gmxpreprocess/ter_db.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/topology/symtab.h"
54 #include "gromacs/utility/cstringutil.h"
55 #include "gromacs/utility/fatalerror.h"
56 #include "gromacs/utility/futil.h"
57 #include "gromacs/utility/smalloc.h"
59 static void copy_atom(t_atoms *atoms1, int a1, t_atoms *atoms2, int a2)
61 atoms2->atom[a2] = atoms1->atom[a1];
62 snew(atoms2->atomname[a2], 1);
63 *atoms2->atomname[a2] = gmx_strdup(*atoms1->atomname[a1]);
66 static int pdbasearch_atom(const char *name, int resind, t_atoms *pdba,
67 const char *searchtype, bool bAllowMissing)
71 for (i = 0; (i < pdba->nr) && (pdba->atom[i].resind != resind); i++)
76 return search_atom(name, i, pdba,
77 searchtype, bAllowMissing);
80 static void hacksearch_atom(int *ii, int *jj, char *name,
81 const int nab[], t_hack *ab[],
82 int resind, t_atoms *pdba)
92 for (i = 0; (i < pdba->nr) && (pdba->atom[i].resind != resind); i++)
96 for (; (i < pdba->nr) && (pdba->atom[i].resind == resind) && (*ii < 0); i++)
98 for (j = 0; (j < nab[i]) && (*ii < 0); j++)
100 if (ab[i][j].nname && strcmp(name, ab[i][j].nname) == 0)
110 static t_hackblock *get_hackblocks(t_atoms *pdba, int nah, t_hackblock ah[],
112 t_hackblock **ntdb, t_hackblock **ctdb,
113 const int *rN, const int *rC)
116 t_hackblock *hb, *ahptr;
119 snew(hb, pdba->nres);
120 /* first the termini */
121 for (i = 0; i < nterpairs; i++)
123 if (ntdb[i] != nullptr)
125 copy_t_hackblock(ntdb[i], &hb[rN[i]]);
127 if (ctdb[i] != nullptr)
129 merge_t_hackblock(ctdb[i], &hb[rC[i]]);
132 /* then the whole hdb */
133 for (rnr = 0; rnr < pdba->nres; rnr++)
135 ahptr = search_h_db(nah, ah, *pdba->resinfo[rnr].rtp);
138 if (hb[rnr].name == nullptr)
140 hb[rnr].name = gmx_strdup(ahptr->name);
142 merge_hacks(ahptr, &hb[rnr]);
148 static void expand_hackblocks_one(t_hackblock *hbr, char *atomname,
149 int *nabi, t_hack **abi, bool bN, bool bC)
154 /* we'll recursively add atoms to atoms */
155 for (j = 0; j < hbr->nhack; j++)
157 /* first check if we're in the N- or C-terminus, then we should ignore
158 all hacks involving atoms from resp. previous or next residue
159 (i.e. which name begins with '-' (N) or '+' (C) */
161 if (bN) /* N-terminus: ignore '-' */
163 for (k = 0; k < 4 && hbr->hack[j].a[k] && !bIgnore; k++)
165 bIgnore = hbr->hack[j].a[k][0] == '-';
168 if (bC) /* C-terminus: ignore '+' */
170 for (k = 0; k < 4 && hbr->hack[j].a[k] && !bIgnore; k++)
172 bIgnore = hbr->hack[j].a[k][0] == '+';
175 /* must be either hdb entry (tp>0) or add from tdb (oname==NULL)
176 and first control aton (AI) matches this atom or
177 delete/replace from tdb (oname!=NULL) and oname matches this atom */
180 ( ( ( hbr->hack[j].tp > 0 || hbr->hack[j].oname == nullptr ) &&
181 strcmp(atomname, hbr->hack[j].ai()) == 0 ) ||
182 ( hbr->hack[j].oname != nullptr &&
183 strcmp(atomname, hbr->hack[j].oname) == 0) ) )
185 /* now expand all hacks for this atom */
186 srenew(*abi, *nabi + hbr->hack[j].nr);
187 for (k = 0; k < hbr->hack[j].nr; k++)
189 copy_t_hack(&hbr->hack[j], &(*abi)[*nabi + k]);
190 (*abi)[*nabi + k].bXSet = FALSE;
191 /* if we're adding (oname==NULL) and don't have a new name (nname)
192 yet, build it from atomname */
193 if ( (*abi)[*nabi + k].nname == nullptr)
195 if ( (*abi)[*nabi + k].oname == nullptr)
197 (*abi)[*nabi + k].nname = gmx_strdup(atomname);
198 (*abi)[*nabi + k].nname[0] = 'H';
205 fprintf(debug, "Hack '%s' %d, replacing nname '%s' with '%s' (old name '%s')\n",
207 (*abi)[*nabi + k].nname, hbr->hack[j].nname,
208 (*abi)[*nabi + k].oname ? (*abi)[*nabi + k].oname : "");
210 sfree((*abi)[*nabi + k].nname);
211 (*abi)[*nabi + k].nname = gmx_strdup(hbr->hack[j].nname);
214 if (hbr->hack[j].tp == 10 && k == 2)
216 /* This is a water virtual site, not a hydrogen */
217 /* Ugly hardcoded name hack */
218 (*abi)[*nabi + k].nname[0] = 'M';
220 else if (hbr->hack[j].tp == 11 && k >= 2)
222 /* This is a water lone pair, not a hydrogen */
223 /* Ugly hardcoded name hack */
224 srenew((*abi)[*nabi + k].nname, 4);
225 (*abi)[*nabi + k].nname[0] = 'L';
226 (*abi)[*nabi + k].nname[1] = 'P';
227 (*abi)[*nabi + k].nname[2] = '1' + k - 2;
228 (*abi)[*nabi + k].nname[3] = '\0';
230 else if (hbr->hack[j].nr > 1)
232 /* adding more than one atom, number them */
233 l = strlen((*abi)[*nabi + k].nname);
234 srenew((*abi)[*nabi + k].nname, l+2);
235 (*abi)[*nabi + k].nname[l] = '1' + k;
236 (*abi)[*nabi + k].nname[l+1] = '\0';
239 (*nabi) += hbr->hack[j].nr;
241 /* add hacks to atoms we've just added */
242 if (hbr->hack[j].tp > 0 || hbr->hack[j].oname == nullptr)
244 for (k = 0; k < hbr->hack[j].nr; k++)
246 expand_hackblocks_one(hbr, (*abi)[*nabi-hbr->hack[j].nr+k].nname,
254 static void expand_hackblocks(t_atoms *pdba, t_hackblock hb[],
255 int nab[], t_hack *ab[],
256 int nterpairs, const int *rN, const int *rC)
261 for (i = 0; i < pdba->nr; i++)
264 for (j = 0; j < nterpairs && !bN; j++)
266 bN = pdba->atom[i].resind == rN[j];
269 for (j = 0; j < nterpairs && !bC; j++)
271 bC = pdba->atom[i].resind == rC[j];
274 /* add hacks to this atom */
275 expand_hackblocks_one(&hb[pdba->atom[i].resind], *pdba->atomname[i],
276 &nab[i], &ab[i], bN, bC);
280 static int check_atoms_present(t_atoms *pdba, const int nab[], t_hack *ab[])
282 int i, j, k, rnr, nadd;
285 for (i = 0; i < pdba->nr; i++)
287 rnr = pdba->atom[i].resind;
288 for (j = 0; j < nab[i]; j++)
290 if (ab[i][j].oname == nullptr)
293 if (ab[i][j].nname == nullptr)
295 gmx_incons("ab[i][j].nname not allocated");
297 /* check if the atom is already present */
298 k = pdbasearch_atom(ab[i][j].nname, rnr, pdba, "check", TRUE);
301 /* We found the added atom. */
302 ab[i][j].bAlreadyPresent = TRUE;
306 ab[i][j].bAlreadyPresent = FALSE;
307 /* count how many atoms we'll add */
311 else if (ab[i][j].nname == nullptr)
322 static void calc_all_pos(t_atoms *pdba, rvec x[], int nab[], t_hack *ab[],
325 int i, j, ii, jj, m, ia, d, rnr, l = 0;
327 rvec xa[4]; /* control atoms for calc_h_pos */
328 rvec xh[MAXH]; /* hydrogen positions from calc_h_pos */
333 for (i = 0; i < pdba->nr; i++)
335 rnr = pdba->atom[i].resind;
336 for (j = 0; j < nab[i]; j += ab[i][j].nr)
338 /* check if we're adding: */
339 if (ab[i][j].oname == nullptr && ab[i][j].tp > 0)
342 for (m = 0; (m < ab[i][j].nctl && bFoundAll); m++)
344 ia = pdbasearch_atom(ab[i][j].a[m], rnr, pdba,
345 bCheckMissing ? "atom" : "check",
349 /* not found in original atoms, might still be in t_hack (ab) */
350 hacksearch_atom(&ii, &jj, ab[i][j].a[m], nab, ab, rnr, pdba);
353 copy_rvec(ab[ii][jj].newx, xa[m]);
360 gmx_fatal(FARGS, "Atom %s not found in residue %s %d"
362 " while adding hydrogens",
364 *pdba->resinfo[rnr].name,
365 pdba->resinfo[rnr].nr,
366 *pdba->resinfo[rnr].rtp);
372 copy_rvec(x[ia], xa[m]);
377 for (m = 0; (m < MAXH); m++)
379 for (d = 0; d < DIM; d++)
391 calc_h_pos(ab[i][j].tp, xa, xh, &l);
392 for (m = 0; m < ab[i][j].nr; m++)
394 copy_rvec(xh[m], ab[i][j+m].newx);
395 ab[i][j+m].bXSet = TRUE;
403 static void free_ab(int natoms, int *nab, t_hack **ab)
407 for (i = 0; i < natoms; i++)
409 free_t_hack(nab[i], &ab[i]);
415 static int add_h_low(t_atoms **pdbaptr, rvec *xptr[],
416 int nah, t_hackblock ah[],
417 int nterpairs, t_hackblock **ntdb, t_hackblock **ctdb,
418 int *rN, int *rC, bool bCheckMissing,
419 int **nabptr, t_hack ***abptr,
420 bool bUpdate_pdba, bool bKeep_old_pdba)
422 t_atoms *newpdba = nullptr, *pdba = nullptr;
424 int i, newi, j, natoms, nalreadypresent;
426 t_hack **ab = nullptr;
431 /* set flags for adding hydrogens (according to hdb) */
437 /* the first time these will be pointers to NULL, but we will
438 return in them the completed arrays, which we will get back
451 /* WOW, everything was already figured out */
452 bUpdate_pdba = FALSE;
456 /* We'll have to do all the hard work */
458 /* first get all the hackblocks for each residue: */
459 hb = get_hackblocks(pdba, nah, ah, nterpairs, ntdb, ctdb, rN, rC);
461 /* expand the hackblocks to atom level */
464 expand_hackblocks(pdba, hb, nab, ab, nterpairs, rN, rC);
465 free_t_hackblock(pdba->nres, &hb);
468 /* Now calc the positions */
469 calc_all_pos(pdba, *xptr, nab, ab, bCheckMissing);
473 /* we don't have to add atoms that are already present in pdba,
474 so we will remove them from the ab (t_hack) */
475 nadd = check_atoms_present(pdba, nab, ab);
477 /* Copy old atoms, making space for new ones */
479 init_t_atoms(newpdba, natoms+nadd, FALSE);
480 newpdba->nres = pdba->nres;
481 sfree(newpdba->resinfo);
482 newpdba->resinfo = pdba->resinfo;
491 /* There is nothing to do: return now */
494 free_ab(natoms, nab, ab);
500 snew(xn, natoms+nadd);
502 for (i = 0; (i < natoms); i++)
504 /* check if this atom wasn't scheduled for deletion */
505 if (nab[i] == 0 || (ab[i][0].nname != nullptr) )
507 if (newi >= natoms+nadd)
509 /*gmx_fatal(FARGS,"Not enough space for adding atoms");*/
511 srenew(xn, natoms+nadd);
514 srenew(newpdba->atom, natoms+nadd);
515 srenew(newpdba->atomname, natoms+nadd);
520 copy_atom(pdba, i, newpdba, newi);
522 copy_rvec((*xptr)[i], xn[newi]);
523 /* process the hacks for this atom */
525 for (j = 0; j < nab[i]; j++)
527 if (ab[i][j].oname == nullptr) /* add */
530 if (newi >= natoms+nadd)
532 /* gmx_fatal(FARGS,"Not enough space for adding atoms");*/
534 srenew(xn, natoms+nadd);
537 srenew(newpdba->atom, natoms+nadd);
538 srenew(newpdba->atomname, natoms+nadd);
543 newpdba->atom[newi].resind = pdba->atom[i].resind;
546 if (ab[i][j].nname != nullptr &&
547 (ab[i][j].oname == nullptr ||
548 strcmp(ab[i][j].oname, *newpdba->atomname[newi]) == 0))
551 if (ab[i][j].oname == nullptr && ab[i][j].bAlreadyPresent)
553 /* This atom is already present, copy it from the input. */
557 copy_atom(pdba, i+nalreadypresent, newpdba, newi);
559 copy_rvec((*xptr)[i+nalreadypresent], xn[newi]);
567 fprintf(debug, "Replacing %d '%s' with (old name '%s') %s\n",
569 (newpdba->atomname[newi] && *newpdba->atomname[newi]) ? *newpdba->atomname[newi] : "",
570 ab[i][j].oname ? ab[i][j].oname : "",
573 snew(newpdba->atomname[newi], 1);
574 *newpdba->atomname[newi] = gmx_strdup(ab[i][j].nname);
575 if (ab[i][j].oname != nullptr && ab[i][j].atom) /* replace */
576 { /* newpdba->atom[newi].m = ab[i][j].atom->m; */
577 /* newpdba->atom[newi].q = ab[i][j].atom->q; */
578 /* newpdba->atom[newi].type = ab[i][j].atom->type; */
583 copy_rvec(ab[i][j].newx, xn[newi]);
586 if (bUpdate_pdba && debug)
588 fprintf(debug, " %s %g %g", *newpdba->atomname[newi],
589 newpdba->atom[newi].m, newpdba->atom[newi].q);
594 i += nalreadypresent;
610 free_ab(natoms, nab, ab);
617 for (i = 0; i < natoms; i++)
619 /* Do not free the atomname string itself, it might be in symtab */
620 /* sfree(*(pdba->atomname[i])); */
621 /* sfree(pdba->atomname[i]); */
623 sfree(pdba->atomname);
625 sfree(pdba->pdbinfo);
637 int add_h(t_atoms **pdbaptr, rvec *xptr[],
638 int nah, t_hackblock ah[],
639 int nterpairs, t_hackblock **ntdb, t_hackblock **ctdb,
640 int *rN, int *rC, bool bAllowMissing,
641 int **nabptr, t_hack ***abptr,
642 bool bUpdate_pdba, bool bKeep_old_pdba)
644 int nold, nnew, niter;
646 /* Here we loop to be able to add atoms to added atoms.
647 * We should not check for missing atoms here.
654 nnew = add_h_low(pdbaptr, xptr, nah, ah, nterpairs, ntdb, ctdb, rN, rC, FALSE,
655 nabptr, abptr, bUpdate_pdba, bKeep_old_pdba);
659 gmx_fatal(FARGS, "More than 100 iterations of add_h. Maybe you are trying to replace an added atom (this is not supported)?");
666 /* Call add_h_low once more, now only for the missing atoms check */
667 add_h_low(pdbaptr, xptr, nah, ah, nterpairs, ntdb, ctdb, rN, rC, TRUE,
668 nabptr, abptr, bUpdate_pdba, bKeep_old_pdba);