src/gromacs/gmxpreprocess/toputil.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) 2012,2014,2015,2017,2018,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.
  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 #include "gmxpre.h"
  38
  39 #include "toputil.h"
  40
  41 #include <climits>
  42 #include <cmath>
  43 #include <cstring>
  44
  45 #include <algorithm>
  46
  47 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
  48 #include "gromacs/gmxpreprocess/notset.h"
  49 #include "gromacs/gmxpreprocess/topdirs.h"
  50 #include "gromacs/topology/block.h"
  51 #include "gromacs/topology/ifunc.h"
  52 #include "gromacs/topology/symtab.h"
  53 #include "gromacs/utility/fatalerror.h"
  54 #include "gromacs/utility/gmxassert.h"
  55 #include "gromacs/utility/smalloc.h"
  56
  57 /* UTILITIES */
  58
  59 void set_p_string(t_param *p, const char *s)
  60 {
  61     if (s)
  62     {
  63         if (strlen(s) < sizeof(p->s)-1)
  64         {
  65             strncpy(p->s, s, sizeof(p->s));
  66         }
  67         else
  68         {
  69             gmx_fatal(FARGS, "Increase MAXSLEN in the grompp code to at least %zu,"
  70                       " or shorten your definition of bonds like %s to at most %d",
  71                       strlen(s)+1, s, MAXSLEN-1);
  72         }
  73     }
  74     else
  75     {
  76         strcpy(p->s, "");
  77     }
  78 }
  79
  80 void pr_alloc (int extra, t_params *pr)
  81 {
  82     int i, j;
  83
  84     /* get new space for arrays */
  85     if (extra < 0)
  86     {
  87         gmx_fatal(FARGS, "Trying to make array smaller.\n");
  88     }
  89     if (extra == 0)
  90     {
  91         return;
  92     }
  93     GMX_ASSERT(pr->nr != 0 || pr->param == nullptr, "Invalid t_params object");
  94     if (pr->nr+extra > pr->maxnr)
  95     {
  96         pr->maxnr = std::max(static_cast<int>(1.2*pr->maxnr), pr->maxnr + extra);
  97         srenew(pr->param, pr->maxnr);
  98         for (i = pr->nr; (i < pr->maxnr); i++)
  99         {
 100             for (j = 0; (j < MAXATOMLIST); j++)
 101             {
 102                 pr->param[i].a[j] = 0;
 103             }
 104             for (j = 0; (j < MAXFORCEPARAM); j++)
 105             {
 106                 pr->param[i].c[j] = 0;
 107             }
 108             set_p_string(&(pr->param[i]), "");
 109         }
 110     }
 111 }
 112
 113 void init_plist(t_params plist[])
 114 {
 115     int i;
 116
 117     for (i = 0; (i < F_NRE); i++)
 118     {
 119         plist[i].nr    = 0;
 120         plist[i].maxnr = 0;
 121         plist[i].param = nullptr;
 122
 123         /* CMAP */
 124         plist[i].ncmap        = 0;
 125         plist[i].cmap         = nullptr;
 126         plist[i].grid_spacing = 0;
 127         plist[i].nc           = 0;
 128         plist[i].nct          = 0;
 129         plist[i].cmap_types   = nullptr;
 130     }
 131 }
 132
 133 void done_plist(t_params *plist)
 134 {
 135     for (int i = 0; i < F_NRE; i++)
 136     {
 137         t_params *pl = &plist[i];
 138         sfree(pl->param);
 139         sfree(pl->cmap);
 140         sfree(pl->cmap_types);
 141     }
 142 }
 143
 144 void cp_param(t_param *dest, t_param *src)
 145 {
 146     int j;
 147
 148     for (j = 0; (j < MAXATOMLIST); j++)
 149     {
 150         dest->a[j] = src->a[j];
 151     }
 152     for (j = 0; (j < MAXFORCEPARAM); j++)
 153     {
 154         dest->c[j] = src->c[j];
 155     }
 156     strncpy(dest->s, src->s, sizeof(dest->s));
 157 }
 158
 159 void add_param_to_list(t_params *list, t_param *b)
 160 {
 161     int j;
 162
 163     /* allocate one position extra */
 164     pr_alloc (1, list);
 165
 166     /* fill the arrays */
 167     for (j = 0; (j < MAXFORCEPARAM); j++)
 168     {
 169         list->param[list->nr].c[j]   = b->c[j];
 170     }
 171     for (j = 0; (j < MAXATOMLIST); j++)
 172     {
 173         list->param[list->nr].a[j]   = b->a[j];
 174     }
 175     memset(list->param[list->nr].s, 0, sizeof(list->param[list->nr].s));
 176
 177     list->nr++;
 178 }
 179
 180
 181 void init_molinfo(t_molinfo *mol)
 182 {
 183     mol->nrexcl     = 0;
 184     mol->excl_set   = FALSE;
 185     mol->bProcessed = FALSE;
 186     init_plist(mol->plist);
 187     init_block(&mol->cgs);
 188     init_block(&mol->mols);
 189     init_blocka(&mol->excls);
 190     init_t_atoms(&mol->atoms, 0, FALSE);
 191 }
 192
 193 /* FREEING MEMORY */
 194
 195 void done_mi(t_molinfo *mi)
 196 {
 197     done_atom (&(mi->atoms));
 198     done_block(&(mi->cgs));
 199     done_block(&(mi->mols));
 200     done_plist(mi->plist);
 201 }
 202
 203 /* PRINTING STRUCTURES */
 204
 205 static void print_bt(FILE *out, Directive d, gpp_atomtype_t at,
 206                      int ftype, int fsubtype, t_params plist[],
 207                      bool bFullDih)
 208 {
 209     /* This dihp is a DIRTY patch because the dih-types do not use
 210      * all four atoms to determine the type.
 211      */
 212     const int    dihp[2][2] = { { 1, 2 }, { 0, 3 } };
 213     t_params    *bt;
 214     int          i, j, f, nral, nrfp;
 215     bool         bDih = FALSE, bSwapParity;
 216
 217     bt = &(plist[ftype]);
 218
 219     if (!bt->nr)
 220     {
 221         return;
 222     }
 223
 224     f = 0;
 225     switch (ftype)
 226     {
 227         case F_G96ANGLES:
 228             f = 1;
 229             break;
 230         case F_G96BONDS:
 231             f = 1;
 232             break;
 233         case F_MORSE:
 234             f = 2;
 235             break;
 236         case F_CUBICBONDS:
 237             f = 3;
 238             break;
 239         case F_CONNBONDS:
 240             f = 4;
 241             break;
 242         case F_HARMONIC:
 243             f = 5;
 244             break;
 245         case F_CROSS_BOND_ANGLES:
 246             f = 2;
 247             break;
 248         case F_CROSS_BOND_BONDS:
 249             f = 3;
 250             break;
 251         case F_UREY_BRADLEY:
 252             f = 4;
 253             break;
 254         case F_PDIHS:
 255         case F_RBDIHS:
 256         case F_FOURDIHS:
 257             bDih = TRUE;
 258             break;
 259         case F_IDIHS:
 260             f    = 1;
 261             bDih = TRUE;
 262             break;
 263         case F_CONSTRNC:
 264             f = 1;
 265             break;
 266         case F_VSITE3FD:
 267             f = 1;
 268             break;
 269         case F_VSITE3FAD:
 270             f = 2;
 271             break;
 272         case F_VSITE3OUT:
 273             f = 3;
 274             break;
 275         case F_VSITE4FDN:
 276             f = 1;
 277             break;
 278         case F_CMAP:
 279             f = 1;
 280             break;
 281
 282         default:
 283             bDih = FALSE;
 284     }
 285     if (bFullDih)
 286     {
 287         bDih = FALSE;
 288     }
 289     if (fsubtype)
 290     {
 291         f = fsubtype-1;
 292     }
 293
 294     nral = NRAL(ftype);
 295     nrfp = NRFP(ftype);
 296
 297     /* header */
 298     fprintf(out, "[ %s ]\n", dir2str(d));
 299     fprintf(out, "; ");
 300     if (!bDih)
 301     {
 302         fprintf (out, "%3s  %4s", "ai", "aj");
 303         for (j = 2; (j < nral); j++)
 304         {
 305             fprintf (out, "  %3c%c", 'a', 'i'+j);
 306         }
 307     }
 308     else
 309     {
 310         for (j = 0; (j < 2); j++)
 311         {
 312             fprintf (out, "%3c%c", 'a', 'i'+dihp[f][j]);
 313         }
 314     }
 315
 316     fprintf (out, " funct");
 317     for (j = 0; (j < nrfp); j++)
 318     {
 319         fprintf (out, " %12c%1d", 'c', j);
 320     }
 321     fprintf (out, "\n");
 322
 323     /* print bondtypes */
 324     for (i = 0; (i < bt->nr); i++)
 325     {
 326         bSwapParity = (bt->param[i].c0() == NOTSET) && (bt->param[i].c1() == -1);
 327         if (!bDih)
 328         {
 329             for (j = 0; (j < nral); j++)
 330             {
 331                 fprintf (out, "%5s ", get_atomtype_name(bt->param[i].a[j], at));
 332             }
 333         }
 334         else
 335         {
 336             for (j = 0; (j < 2); j++)
 337             {
 338                 fprintf (out, "%5s ", get_atomtype_name(bt->param[i].a[dihp[f][j]], at));
 339             }
 340         }
 341         fprintf (out, "%5d ", bSwapParity ? -f-1 : f+1);
 342
 343         if (bt->param[i].s[0])
 344         {
 345             fprintf(out, "   %s", bt->param[i].s);
 346         }
 347         else
 348         {
 349             for (j = 0; (j < nrfp && (bt->param[i].c[j] != NOTSET)); j++)
 350             {
 351                 fprintf (out, "%13.6e ", bt->param[i].c[j]);
 352             }
 353         }
 354
 355         fprintf (out, "\n");
 356     }
 357     fprintf (out, "\n");
 358     fflush (out);
 359 }
 360
 361 void print_blocka(FILE *out, const char *szName,
 362                   const char *szIndex, const char *szA,
 363                   t_blocka *block)
 364 {
 365     int i, j;
 366
 367     fprintf (out, "; %s\n", szName);
 368     fprintf (out, "; %4s    %s\n", szIndex, szA);
 369     for (i = 0; (i < block->nr); i++)
 370     {
 371         for (i = 0; (i < block->nr); i++)
 372         {
 373             fprintf (out, "%6d", i+1);
 374             for (j = block->index[i]; (j < (block->index[i+1])); j++)
 375             {
 376                 fprintf (out, "%5d", block->a[j]+1);
 377             }
 378             fprintf (out, "\n");
 379         }
 380         fprintf (out, "\n");
 381     }
 382 }
 383
 384 void print_excl(FILE *out, int natoms, t_excls excls[])
 385 {
 386     int         i;
 387     bool        have_excl;
 388     int         j;
 389
 390     have_excl = FALSE;
 391     for (i = 0; i < natoms && !have_excl; i++)
 392     {
 393         have_excl = (excls[i].nr > 0);
 394     }
 395
 396     if (have_excl)
 397     {
 398         fprintf (out, "[ %s ]\n", dir2str(Directive::d_exclusions));
 399         fprintf (out, "; %4s    %s\n", "i", "excluded from i");
 400         for (i = 0; i < natoms; i++)
 401         {
 402             if (excls[i].nr > 0)
 403             {
 404                 fprintf (out, "%6d ", i+1);
 405                 for (j = 0; j < excls[i].nr; j++)
 406                 {
 407                     fprintf (out, " %5d", excls[i].e[j]+1);
 408                 }
 409                 fprintf (out, "\n");
 410             }
 411         }
 412         fprintf (out, "\n");
 413         fflush(out);
 414     }
 415 }
 416
 417 static double get_residue_charge(const t_atoms *atoms, int at)
 418 {
 419     int    ri;
 420     double q;
 421
 422     ri = atoms->atom[at].resind;
 423     q  = 0;
 424     while (at < atoms->nr && atoms->atom[at].resind == ri)
 425     {
 426         q += atoms->atom[at].q;
 427         at++;
 428     }
 429
 430     return q;
 431 }
 432
 433 void print_atoms(FILE *out, gpp_atomtype_t atype, t_atoms *at, int *cgnr,
 434                  bool bRTPresname)
 435 {
 436     int         i, ri;
 437     int         tpA, tpB;
 438     const char *as;
 439     char       *tpnmA, *tpnmB;
 440     double      qres, qtot;
 441
 442     as = dir2str(Directive::d_atoms);
 443     fprintf(out, "[ %s ]\n", as);
 444     fprintf(out, "; %4s %10s %6s %7s%6s %6s %10s %10s %6s %10s %10s\n",
 445             "nr", "type", "resnr", "residue", "atom", "cgnr", "charge", "mass", "typeB", "chargeB", "massB");
 446
 447     qtot  = 0;
 448
 449     if (at->nres)
 450     {
 451         /* if the information is present... */
 452         for (i = 0; (i < at->nr); i++)
 453         {
 454             ri = at->atom[i].resind;
 455             if ((i == 0 || ri != at->atom[i-1].resind) &&
 456                 at->resinfo[ri].rtp != nullptr)
 457             {
 458                 qres = get_residue_charge(at, i);
 459                 fprintf(out, "; residue %3d %-3s rtp %-4s q ",
 460                         at->resinfo[ri].nr,
 461                         *at->resinfo[ri].name,
 462                         *at->resinfo[ri].rtp);
 463                 if (fabs(qres) < 0.001)
 464                 {
 465                     fprintf(out, " %s", "0.0");
 466                 }
 467                 else
 468                 {
 469                     fprintf(out, "%+3.1f", qres);
 470                 }
 471                 fprintf(out, "\n");
 472             }
 473             tpA = at->atom[i].type;
 474             if ((tpnmA = get_atomtype_name(tpA, atype)) == nullptr)
 475             {
 476                 gmx_fatal(FARGS, "tpA = %d, i= %d in print_atoms", tpA, i);
 477             }
 478
 479             /* This is true by construction, but static analysers don't know */
 480             GMX_ASSERT(!bRTPresname || at->resinfo[at->atom[i].resind].rtp, "-rtpres did not have residue name available");
 481             fprintf(out, "%6d %10s %6d%c %5s %6s %6d %10g %10g",
 482                     i+1, tpnmA,
 483                     at->resinfo[ri].nr,
 484                     at->resinfo[ri].ic,
 485                     bRTPresname ?
 486                     *(at->resinfo[at->atom[i].resind].rtp) :
 487                     *(at->resinfo[at->atom[i].resind].name),
 488                     *(at->atomname[i]), cgnr[i],
 489                     at->atom[i].q, at->atom[i].m);
 490             if (PERTURBED(at->atom[i]))
 491             {
 492                 tpB = at->atom[i].typeB;
 493                 if ((tpnmB = get_atomtype_name(tpB, atype)) == nullptr)
 494                 {
 495                     gmx_fatal(FARGS, "tpB = %d, i= %d in print_atoms", tpB, i);
 496                 }
 497                 fprintf(out, " %6s %10g %10g",
 498                         tpnmB, at->atom[i].qB, at->atom[i].mB);
 499             }
 500             // Accumulate the total charge to help troubleshoot issues.
 501             qtot += static_cast<double>(at->atom[i].q);
 502             // Round it to zero if it is close to zero, because
 503             // printing -9.34e-5 confuses users.
 504             if (fabs(qtot) < 0.0001)
 505             {
 506                 qtot = 0;
 507             }
 508             // Write the total charge for the last atom of the system
 509             // and/or residue, because generally that's where it is
 510             // expected to be an integer.
 511             if (i == at->nr-1 || ri != at->atom[i+1].resind)
 512             {
 513                 fprintf(out, "   ; qtot %.4g\n", qtot);
 514             }
 515             else
 516             {
 517                 fputs("\n", out);
 518             }
 519         }
 520     }
 521     fprintf(out, "\n");
 522     fflush(out);
 523 }
 524
 525 void print_bondeds(FILE *out, int natoms, Directive d,
 526                    int ftype, int fsubtype, t_params plist[])
 527 {
 528     t_symtab       stab;
 529     gpp_atomtype_t atype;
 530     t_param       *param;
 531     t_atom        *a;
 532     int            i;
 533
 534     atype = init_atomtype();
 535     snew(a, 1);
 536     snew(param, 1);
 537     open_symtab(&stab);
 538     for (i = 0; (i < natoms); i++)
 539     {
 540         char buf[12];
 541         sprintf(buf, "%4d", (i+1));
 542         add_atomtype(atype, &stab, a, buf, param, 0, 0);
 543     }
 544     print_bt(out, d, atype, ftype, fsubtype, plist, TRUE);
 545
 546     done_symtab(&stab);
 547     sfree(a);
 548     sfree(param);
 549     done_atomtype(atype);
 550 }