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