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 add_param_to_list(InteractionTypeParameters *list, const InteractionType &b)
  61 {
  62     list->interactionTypes.emplace_back(b);
  63 }
  64
  65 /* PRINTING STRUCTURES */
  66
  67 static void print_bt(FILE *out, Directive d, PreprocessingAtomTypes *at,
  68                      int ftype, int fsubtype, gmx::ArrayRef<const InteractionTypeParameters> plist,
  69                      bool bFullDih)
  70 {
  71     /* This dihp is a DIRTY patch because the dih-types do not use
  72      * all four atoms to determine the type.
  73      */
  74     const int                        dihp[2][2] = { { 1, 2 }, { 0, 3 } };
  75     int                              nral, nrfp;
  76     bool                             bDih = false, bSwapParity;
  77
  78     const InteractionTypeParameters *bt = &(plist[ftype]);
  79
  80     if (bt->size() == 0)
  81     {
  82         return;
  83     }
  84
  85     int f = 0;
  86     switch (ftype)
  87     {
  88         case F_G96ANGLES:
  89             f = 1;
  90             break;
  91         case F_G96BONDS:
  92             f = 1;
  93             break;
  94         case F_MORSE:
  95             f = 2;
  96             break;
  97         case F_CUBICBONDS:
  98             f = 3;
  99             break;
 100         case F_CONNBONDS:
 101             f = 4;
 102             break;
 103         case F_HARMONIC:
 104             f = 5;
 105             break;
 106         case F_CROSS_BOND_ANGLES:
 107             f = 2;
 108             break;
 109         case F_CROSS_BOND_BONDS:
 110             f = 3;
 111             break;
 112         case F_UREY_BRADLEY:
 113             f = 4;
 114             break;
 115         case F_PDIHS:
 116         case F_RBDIHS:
 117         case F_FOURDIHS:
 118             bDih = TRUE;
 119             break;
 120         case F_IDIHS:
 121             f    = 1;
 122             bDih = TRUE;
 123             break;
 124         case F_CONSTRNC:
 125             f = 1;
 126             break;
 127         case F_VSITE3FD:
 128             f = 1;
 129             break;
 130         case F_VSITE3FAD:
 131             f = 2;
 132             break;
 133         case F_VSITE3OUT:
 134             f = 3;
 135             break;
 136         case F_VSITE4FDN:
 137             f = 1;
 138             break;
 139         case F_CMAP:
 140             f = 1;
 141             break;
 142
 143         default:
 144             bDih = FALSE;
 145     }
 146     if (bFullDih)
 147     {
 148         bDih = FALSE;
 149     }
 150     if (fsubtype)
 151     {
 152         f = fsubtype-1;
 153     }
 154
 155     nral = NRAL(ftype);
 156     nrfp = NRFP(ftype);
 157
 158     /* header */
 159     fprintf(out, "[ %s ]\n", dir2str(d));
 160     fprintf(out, "; ");
 161     if (!bDih)
 162     {
 163         fprintf (out, "%3s  %4s", "ai", "aj");
 164         for (int j = 2; (j < nral); j++)
 165         {
 166             fprintf (out, "  %3c%c", 'a', 'i'+j);
 167         }
 168     }
 169     else
 170     {
 171         for (int j = 0; (j < 2); j++)
 172         {
 173             fprintf (out, "%3c%c", 'a', 'i'+dihp[f][j]);
 174         }
 175     }
 176
 177     fprintf (out, " funct");
 178     for (int j = 0; (j < nrfp); j++)
 179     {
 180         fprintf (out, " %12c%1d", 'c', j);
 181     }
 182     fprintf (out, "\n");
 183
 184     /* print bondtypes */
 185     for (const auto &parm : bt->interactionTypes)
 186     {
 187         bSwapParity = (parm.c0() == NOTSET) && (parm.c1() == -1);
 188         gmx::ArrayRef<const int> atoms = parm.atoms();
 189         if (!bDih)
 190         {
 191             for (int j = 0; (j < nral); j++)
 192             {
 193                 fprintf (out, "%5s ", at->atomNameFromAtomType(atoms[j]));
 194             }
 195         }
 196         else
 197         {
 198             for (int j = 0; (j < 2); j++)
 199             {
 200                 fprintf (out, "%5s ", at->atomNameFromAtomType(atoms[dihp[f][j]]));
 201             }
 202         }
 203         fprintf (out, "%5d ", bSwapParity ? -f-1 : f+1);
 204
 205         if (!parm.interactionTypeName().empty())
 206         {
 207             fprintf(out, "   %s", parm.interactionTypeName().c_str());
 208         }
 209         else
 210         {
 211             gmx::ArrayRef<const real> forceParam = parm.forceParam();
 212             for (int j = 0; (j < nrfp) && (forceParam[j] != NOTSET); j++)
 213             {
 214                 fprintf (out, "%13.6e ", forceParam[j]);
 215             }
 216         }
 217
 218         fprintf (out, "\n");
 219     }
 220     fprintf (out, "\n");
 221     fflush (out);
 222 }
 223
 224 void print_blocka(FILE *out, const char *szName,
 225                   const char *szIndex, const char *szA,
 226                   t_blocka *block)
 227 {
 228     int i, j;
 229
 230     fprintf (out, "; %s\n", szName);
 231     fprintf (out, "; %4s    %s\n", szIndex, szA);
 232     for (i = 0; (i < block->nr); i++)
 233     {
 234         for (i = 0; (i < block->nr); i++)
 235         {
 236             fprintf (out, "%6d", i+1);
 237             for (j = block->index[i]; (j < (block->index[i+1])); j++)
 238             {
 239                 fprintf (out, "%5d", block->a[j]+1);
 240             }
 241             fprintf (out, "\n");
 242         }
 243         fprintf (out, "\n");
 244     }
 245 }
 246
 247 void print_excl(FILE *out, int natoms, t_excls excls[])
 248 {
 249     int         i;
 250     bool        have_excl;
 251     int         j;
 252
 253     have_excl = FALSE;
 254     for (i = 0; i < natoms && !have_excl; i++)
 255     {
 256         have_excl = (excls[i].nr > 0);
 257     }
 258
 259     if (have_excl)
 260     {
 261         fprintf (out, "[ %s ]\n", dir2str(Directive::d_exclusions));
 262         fprintf (out, "; %4s    %s\n", "i", "excluded from i");
 263         for (i = 0; i < natoms; i++)
 264         {
 265             if (excls[i].nr > 0)
 266             {
 267                 fprintf (out, "%6d ", i+1);
 268                 for (j = 0; j < excls[i].nr; j++)
 269                 {
 270                     fprintf (out, " %5d", excls[i].e[j]+1);
 271                 }
 272                 fprintf (out, "\n");
 273             }
 274         }
 275         fprintf (out, "\n");
 276         fflush(out);
 277     }
 278 }
 279
 280 static double get_residue_charge(const t_atoms *atoms, int at)
 281 {
 282     int    ri;
 283     double q;
 284
 285     ri = atoms->atom[at].resind;
 286     q  = 0;
 287     while (at < atoms->nr && atoms->atom[at].resind == ri)
 288     {
 289         q += atoms->atom[at].q;
 290         at++;
 291     }
 292
 293     return q;
 294 }
 295
 296 void print_atoms(FILE *out, PreprocessingAtomTypes *atype, t_atoms *at, int *cgnr,
 297                  bool bRTPresname)
 298 {
 299     int               i, ri;
 300     int               tpA, tpB;
 301     const char       *as;
 302     const char       *tpnmA, *tpnmB;
 303     double            qres, qtot;
 304
 305     as = dir2str(Directive::d_atoms);
 306     fprintf(out, "[ %s ]\n", as);
 307     fprintf(out, "; %4s %10s %6s %7s%6s %6s %10s %10s %6s %10s %10s\n",
 308             "nr", "type", "resnr", "residue", "atom", "cgnr", "charge", "mass", "typeB", "chargeB", "massB");
 309
 310     qtot  = 0;
 311
 312     if (at->nres)
 313     {
 314         /* if the information is present... */
 315         for (i = 0; (i < at->nr); i++)
 316         {
 317             ri = at->atom[i].resind;
 318             if ((i == 0 || ri != at->atom[i-1].resind) &&
 319                 at->resinfo[ri].rtp != nullptr)
 320             {
 321                 qres = get_residue_charge(at, i);
 322                 fprintf(out, "; residue %3d %-3s rtp %-4s q ",
 323                         at->resinfo[ri].nr,
 324                         *at->resinfo[ri].name,
 325                         *at->resinfo[ri].rtp);
 326                 if (fabs(qres) < 0.001)
 327                 {
 328                     fprintf(out, " %s", "0.0");
 329                 }
 330                 else
 331                 {
 332                     fprintf(out, "%+3.1f", qres);
 333                 }
 334                 fprintf(out, "\n");
 335             }
 336             tpA = at->atom[i].type;
 337             if ((tpnmA = atype->atomNameFromAtomType(tpA)) == nullptr)
 338             {
 339                 gmx_fatal(FARGS, "tpA = %d, i= %d in print_atoms", tpA, i);
 340             }
 341
 342             /* This is true by construction, but static analysers don't know */
 343             GMX_ASSERT(!bRTPresname || at->resinfo[at->atom[i].resind].rtp, "-rtpres did not have residue name available");
 344             fprintf(out, "%6d %10s %6d%c %5s %6s %6d %10g %10g",
 345                     i+1, tpnmA,
 346                     at->resinfo[ri].nr,
 347                     at->resinfo[ri].ic,
 348                     bRTPresname ?
 349                     *(at->resinfo[at->atom[i].resind].rtp) :
 350                     *(at->resinfo[at->atom[i].resind].name),
 351                     *(at->atomname[i]), cgnr[i],
 352                     at->atom[i].q, at->atom[i].m);
 353             if (PERTURBED(at->atom[i]))
 354             {
 355                 tpB = at->atom[i].typeB;
 356                 if ((tpnmB = atype->atomNameFromAtomType(tpB)) == nullptr)
 357                 {
 358                     gmx_fatal(FARGS, "tpB = %d, i= %d in print_atoms", tpB, i);
 359                 }
 360                 fprintf(out, " %6s %10g %10g",
 361                         tpnmB, at->atom[i].qB, at->atom[i].mB);
 362             }
 363             // Accumulate the total charge to help troubleshoot issues.
 364             qtot += static_cast<double>(at->atom[i].q);
 365             // Round it to zero if it is close to zero, because
 366             // printing -9.34e-5 confuses users.
 367             if (fabs(qtot) < 0.0001)
 368             {
 369                 qtot = 0;
 370             }
 371             // Write the total charge for the last atom of the system
 372             // and/or residue, because generally that's where it is
 373             // expected to be an integer.
 374             if (i == at->nr-1 || ri != at->atom[i+1].resind)
 375             {
 376                 fprintf(out, "   ; qtot %.4g\n", qtot);
 377             }
 378             else
 379             {
 380                 fputs("\n", out);
 381             }
 382         }
 383     }
 384     fprintf(out, "\n");
 385     fflush(out);
 386 }
 387
 388 void print_bondeds(FILE *out, int natoms, Directive d,
 389                    int ftype, int fsubtype, gmx::ArrayRef<const InteractionTypeParameters> plist)
 390 {
 391     t_symtab               stab;
 392     t_atom                *a;
 393
 394     PreprocessingAtomTypes atype;
 395     snew(a, 1);
 396     open_symtab(&stab);
 397     for (int i = 0; (i < natoms); i++)
 398     {
 399         char buf[12];
 400         sprintf(buf, "%4d", (i+1));
 401         atype.addType(&stab, *a, buf, InteractionType({}, {}), 0, 0);
 402     }
 403     print_bt(out, d, &atype, ftype, fsubtype, plist, TRUE);
 404
 405     done_symtab(&stab);
 406     sfree(a);
 407 }