src/gromacs/gmxpreprocess/hackblock.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2011,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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 /* This file is completely threadsafe - keep it that way! */
  38 #include "gmxpre.h"
  39
  40 #include "hackblock.h"
  41
  42 #include <cstring>
  43
  44 #include "gromacs/gmxpreprocess/notset.h"
  45 #include "gromacs/math/vec.h"
  46 #include "gromacs/mdtypes/md_enums.h"
  47 #include "gromacs/utility/cstringutil.h"
  48 #include "gromacs/utility/fatalerror.h"
  49 #include "gromacs/utility/smalloc.h"
  50
  51 /* these MUST correspond to the enum in hackblock.h */
  52 const char *btsNames[ebtsNR] = { "bonds", "angles", "dihedrals", "impropers", "exclusions", "cmap" };
  53 const int   btsNiatoms[ebtsNR] = { 2,       3,        4,           4,           2,             5 };
  54
  55 static void free_t_bonded(t_rbonded *rb)
  56 {
  57     int i;
  58
  59     for (i = 0; i < MAXATOMLIST; i++)
  60     {
  61         sfree(rb->a[i]);
  62     }
  63     sfree(rb->s);
  64 }
  65
  66 static void free_t_bondeds(t_rbondeds *rbs)
  67 {
  68     int i;
  69
  70     for (i = 0; i < rbs->nb; i++)
  71     {
  72         free_t_bonded(&rbs->b[i]);
  73     }
  74     sfree(rbs->b);
  75     rbs->b  = nullptr;
  76     rbs->nb = 0;
  77 }
  78
  79 void free_t_restp(int nrtp, t_restp **rtp)
  80 {
  81     int i, j;
  82
  83     for (i = 0; i < nrtp; i++)
  84     {
  85         sfree((*rtp)[i].resname);
  86         sfree((*rtp)[i].atom);
  87         for (j = 0; j < (*rtp)[i].natom; j++)
  88         {
  89             sfree(*(*rtp)[i].atomname[j]);
  90             sfree((*rtp)[i].atomname[j]);
  91         }
  92         sfree((*rtp)[i].atomname);
  93         sfree((*rtp)[i].cgnr);
  94         for (j = 0; j < ebtsNR; j++)
  95         {
  96             free_t_bondeds(&(*rtp)[i].rb[j]);
  97         }
  98     }
  99     sfree(*rtp);
 100 }
 101
 102 void free_t_hack(int nh, t_hack **h)
 103 {
 104     int i, j;
 105
 106     for (i = 0; i < nh; i++)
 107     {
 108         sfree((*h)[i].oname);
 109         sfree((*h)[i].nname);
 110         sfree((*h)[i].atom);
 111         for (j = 0; j < 4; j++)
 112         {
 113             sfree((*h)[i].a[j]);
 114         }
 115     }
 116     sfree(*h);
 117     *h = nullptr;
 118 }
 119
 120 void free_t_hackblock(int nhb, t_hackblock **hb)
 121 {
 122     int i, j;
 123
 124     for (i = 0; i < nhb; i++)
 125     {
 126         sfree((*hb)[i].name);
 127         free_t_hack((*hb)[i].nhack, &(*hb)[i].hack);
 128         for (j = 0; j < ebtsNR; j++)
 129         {
 130             free_t_bondeds(&(*hb)[i].rb[j]);
 131         }
 132     }
 133     sfree(*hb);
 134 }
 135
 136 void clear_t_hackblock(t_hackblock *hb)
 137 {
 138     int i;
 139
 140     hb->name    = nullptr;
 141     hb->nhack   = 0;
 142     hb->maxhack = 0;
 143     hb->hack    = nullptr;
 144     for (i = 0; i < ebtsNR; i++)
 145     {
 146         hb->rb[i].nb = 0;
 147         hb->rb[i].b  = nullptr;
 148     }
 149 }
 150
 151 void clear_t_hack(t_hack *hack)
 152 {
 153     int i;
 154
 155     hack->nr    = 0;
 156     hack->oname = nullptr;
 157     hack->nname = nullptr;
 158     hack->atom  = nullptr;
 159     hack->cgnr  = NOTSET;
 160     hack->tp    = 0;
 161     hack->nctl  = 0;
 162     for (i = 0; i < 4; i++)
 163     {
 164         hack->a[i]  = nullptr;
 165     }
 166     for (i = 0; i < DIM; i++)
 167     {
 168         hack->newx[i] = NOTSET;
 169     }
 170 }
 171
 172 #define safe_strdup(str) (((str) != NULL) ? gmx_strdup(str) : NULL)
 173
 174 static void copy_t_rbonded(t_rbonded *s, t_rbonded *d)
 175 {
 176     int i;
 177
 178     for (i = 0; i < MAXATOMLIST; i++)
 179     {
 180         d->a[i] = safe_strdup(s->a[i]);
 181     }
 182     d->s     = safe_strdup(s->s);
 183     d->match = s->match;
 184 }
 185
 186 static bool contains_char(t_rbonded *s, char c)
 187 {
 188     int      i;
 189     bool     bRet;
 190
 191     bRet = FALSE;
 192     for (i = 0; i < MAXATOMLIST; i++)
 193     {
 194         if (s->a[i] && s->a[i][0] == c)
 195         {
 196             bRet = TRUE;
 197         }
 198     }
 199
 200     return bRet;
 201 }
 202
 203 static int
 204 rbonded_find_atoms_in_list(t_rbonded *b, t_rbonded blist[], int nlist, int natoms)
 205 {
 206     int      i, k;
 207     int      foundPos = -1;
 208     bool     atomsMatch;
 209
 210     for (i = 0; i < nlist && foundPos < 0; i++)
 211     {
 212         atomsMatch = TRUE;
 213         for (k = 0; k < natoms && atomsMatch; k++)
 214         {
 215             atomsMatch = atomsMatch && !strcmp(b->a[k], blist[i].a[k]);
 216         }
 217         /* Try reverse if forward match did not work */
 218         if (!atomsMatch)
 219         {
 220             atomsMatch = TRUE;
 221             for (k = 0; k < natoms && atomsMatch; k++)
 222             {
 223                 atomsMatch = atomsMatch && !strcmp(b->a[k], blist[i].a[natoms-1-k]);
 224             }
 225         }
 226         if (atomsMatch)
 227         {
 228             foundPos = i;
 229             /* If all the atoms AND all the parameters match, it is likely that
 230              * the user made a copy-and-paste mistake (since it would be much cheaper
 231              * to just bump the force constant 2x if you really want it twice).
 232              * Since we only have the unparsed string here we can only detect
 233              * EXACT matches (including identical whitespace).
 234              */
 235             if (!strcmp(b->s, blist[i].s))
 236             {
 237                 gmx_warning("Duplicate line found in or between hackblock and rtp entries");
 238             }
 239         }
 240     }
 241     return foundPos;
 242 }
 243
 244 bool merge_t_bondeds(t_rbondeds s[], t_rbondeds d[], bool bMin, bool bPlus)
 245 {
 246     int      i, j;
 247     bool     bBondsRemoved;
 248     int      nbHackblockStart;
 249     int      index;
 250
 251     bBondsRemoved = FALSE;
 252     for (i = 0; i < ebtsNR; i++)
 253     {
 254         if (s[i].nb > 0)
 255         {
 256             /* Record how many bonds we have in the destination when we start.
 257              *
 258              * If an entry is present in the hackblock (destination), we will
 259              * not add the one from the main rtp, since the point is for hackblocks
 260              * to overwrite it. However, if there is no hackblock entry we do
 261              * allow multiple main rtp entries since some forcefield insist on that.
 262              *
 263              * We accomplish this by checking the position we find an entry in,
 264              * rather than merely checking whether it exists at all.
 265              * If that index is larger than the original (hackblock) destination
 266              * size, it was added from the main rtp, and then we will allow more
 267              * such entries. In contrast, if the entry found has a lower index
 268              * it is a hackblock entry meant to override the main rtp, and then
 269              * we don't add the main rtp one.
 270              */
 271             nbHackblockStart = d[i].nb;
 272
 273             /* make space */
 274             srenew(d[i].b, d[i].nb + s[i].nb);
 275             for (j = 0; j < s[i].nb; j++)
 276             {
 277                 /* Check if this bonded string already exists before adding.
 278                  * We are merging from the main RTP to the hackblocks, so this
 279                  * will mean the hackblocks overwrite the man RTP, as intended.
 280                  */
 281                 index = rbonded_find_atoms_in_list(&s[i].b[j], d[i].b, d[i].nb, btsNiatoms[i]);
 282                 /* - If we did not find this interaction at all, the index will be -1,
 283                  *   and then we should definitely add it to the merged hackblock and rtp.
 284                  *
 285                  * Alternatively, if it was found, index will be >=0.
 286                  * - In case this index is lower than the original number of entries,
 287                  *   it is already present as a *hackblock* entry, and those should
 288                  *   always override whatever we have listed in the RTP. Thus, we
 289                  *   should just keep that one and not add anything from the RTP.
 290                  * - Finally, if it was found, but with an index higher than
 291                  *   the original number of entries, it comes from the RTP rather
 292                  *   than hackblock, and then we must have added it ourselves
 293                  *   in a previous iteration. In that case it is a matter of
 294                  *   several entries for the same sequence of atoms, and we allow
 295                  *   that in the RTP. In this case we should simply copy all of
 296                  *   them, including this one.
 297                  */
 298                 if (index < 0 || index >= nbHackblockStart)
 299                 {
 300                     if (!(bMin && contains_char(&s[i].b[j], '-'))
 301                         && !(bPlus && contains_char(&s[i].b[j], '+')))
 302                     {
 303                         copy_t_rbonded(&s[i].b[j], &d[i].b[ d[i].nb ]);
 304                         d[i].nb++;
 305                     }
 306                     else if (i == ebtsBONDS)
 307                     {
 308                         bBondsRemoved = TRUE;
 309                     }
 310                 }
 311                 else
 312                 {
 313                     /* This is the common case where a hackblock entry simply
 314                      * overrides the RTP, so we cannot warn here.
 315                      */
 316                 }
 317             }
 318         }
 319     }
 320     return bBondsRemoved;
 321 }
 322
 323 void copy_t_restp(t_restp *s, t_restp *d)
 324 {
 325     int i;
 326
 327     *d         = *s;
 328     d->resname = safe_strdup(s->resname);
 329     snew(d->atom, s->natom);
 330     for (i = 0; i < s->natom; i++)
 331     {
 332         d->atom[i] = s->atom[i];
 333     }
 334     snew(d->atomname, s->natom);
 335     for (i = 0; i < s->natom; i++)
 336     {
 337         snew(d->atomname[i], 1);
 338         *d->atomname[i] = safe_strdup(*s->atomname[i]);
 339     }
 340     snew(d->cgnr, s->natom);
 341     for (i = 0; i < s->natom; i++)
 342     {
 343         d->cgnr[i] = s->cgnr[i];
 344     }
 345     for (i = 0; i < ebtsNR; i++)
 346     {
 347         d->rb[i].type = s->rb[i].type;
 348         d->rb[i].nb   = 0;
 349         d->rb[i].b    = nullptr;
 350     }
 351     merge_t_bondeds(s->rb, d->rb, FALSE, FALSE);
 352 }
 353
 354 void copy_t_hack(t_hack *s, t_hack *d)
 355 {
 356     int i;
 357
 358     *d       = *s;
 359     d->oname = safe_strdup(s->oname);
 360     d->nname = safe_strdup(s->nname);
 361     if (s->atom)
 362     {
 363         snew(d->atom, 1);
 364         *(d->atom) = *(s->atom);
 365     }
 366     else
 367     {
 368         d->atom = nullptr;
 369     }
 370     for (i = 0; i < 4; i++)
 371     {
 372         d->a[i] = safe_strdup(s->a[i]);
 373     }
 374     copy_rvec(s->newx, d->newx);
 375 }
 376
 377 void merge_hacks_lo(int ns, t_hack *s, int *nd, t_hack **d)
 378 {
 379     int i;
 380
 381     if (ns)
 382     {
 383         srenew(*d, *nd + ns);
 384         for (i = 0; i < ns; i++)
 385         {
 386             copy_t_hack(&s[i], &(*d)[*nd + i]);
 387         }
 388         (*nd) += ns;
 389     }
 390 }
 391
 392 void merge_hacks(t_hackblock *s, t_hackblock *d)
 393 {
 394     merge_hacks_lo(s->nhack, s->hack, &d->nhack, &d->hack);
 395 }
 396
 397 void merge_t_hackblock(t_hackblock *s, t_hackblock *d)
 398 {
 399     merge_hacks(s, d);
 400     merge_t_bondeds(s->rb, d->rb, FALSE, FALSE);
 401 }
 402
 403 void copy_t_hackblock(t_hackblock *s, t_hackblock *d)
 404 {
 405     int i;
 406
 407     *d       = *s;
 408     d->name  = safe_strdup(s->name);
 409     d->nhack = 0;
 410     d->hack  = nullptr;
 411     for (i = 0; i < ebtsNR; i++)
 412     {
 413         d->rb[i].nb = 0;
 414         d->rb[i].b  = nullptr;
 415     }
 416     merge_t_hackblock(s, d);
 417 }
 418
 419 #undef safe_strdup
 420
 421 void dump_hb(FILE *out, int nres, t_hackblock hb[])
 422 {
 423     int i, j, k, l;
 424
 425 #define SS(s) (s) ? (s) : "-"
 426 #define SA(s) (s) ? "+" : ""
 427     fprintf(out, "t_hackblock\n");
 428     for (i = 0; i < nres; i++)
 429     {
 430         fprintf(out, "%3d %4s %2d %2d\n",
 431                 i, SS(hb[i].name), hb[i].nhack, hb[i].maxhack);
 432         if (hb[i].nhack)
 433         {
 434             for (j = 0; j < hb[i].nhack; j++)
 435             {
 436                 fprintf(out, "%d: %d %4s %4s %1s %2d %d %4s %4s %4s %4s\n",
 437                         j, hb[i].hack[j].nr,
 438                         SS(hb[i].hack[j].oname), SS(hb[i].hack[j].nname),
 439                         SA(hb[i].hack[j].atom), hb[i].hack[j].tp, hb[i].hack[j].cgnr,
 440                         SS(hb[i].hack[j].ai()), SS(hb[i].hack[j].aj()),
 441                         SS(hb[i].hack[j].ak()), SS(hb[i].hack[j].al()) );
 442             }
 443         }
 444         for (j = 0; j < ebtsNR; j++)
 445         {
 446             if (hb[i].rb[j].nb)
 447             {
 448                 fprintf(out, " %c %d:", btsNames[j][0], hb[i].rb[j].nb);
 449                 for (k = 0; k < hb[i].rb[j].nb; k++)
 450                 {
 451                     fprintf(out, " [");
 452                     for (l = 0; l < btsNiatoms[j]; l++)
 453                     {
 454                         fprintf(out, " %s", hb[i].rb[j].b[k].a[l]);
 455                     }
 456                     fprintf(out, " %s]", SS(hb[i].rb[j].b[k].s));
 457                 }
 458                 fprintf(out, " Entry matched: %s\n", yesno_names[hb[i].rb[j].b[k].match]);
 459             }
 460         }
 461         fprintf(out, "\n");
 462     }
 463 #undef SS
 464 #undef SA
 465 }