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 by the GROMACS development team.
   7  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "toputil.h"
  41
  42 #include <climits>
  43 #include <cmath>
  44 #include <cstring>
  45
  46 #include <algorithm>
  47
  48 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
  49 #include "gromacs/gmxpreprocess/grompp_impl.h"
  50 #include "gromacs/gmxpreprocess/notset.h"
  51 #include "gromacs/gmxpreprocess/topdirs.h"
  52 #include "gromacs/topology/block.h"
  53 #include "gromacs/topology/ifunc.h"
  54 #include "gromacs/topology/symtab.h"
  55 #include "gromacs/utility/fatalerror.h"
  56 #include "gromacs/utility/gmxassert.h"
  57 #include "gromacs/utility/smalloc.h"
  58
  59 /* UTILITIES */
  60
  61 void add_param_to_list(InteractionsOfType* list, const InteractionOfType& b)
  62 {
  63     list->interactionTypes.emplace_back(b);
  64 }
  65
  66 /* PRINTING STRUCTURES */
  67
  68 static void print_bt(FILE*                                   out,
  69                      Directive                               d,
  70                      PreprocessingAtomTypes*                 at,
  71                      int                                     ftype,
  72                      int                                     fsubtype,
  73                      gmx::ArrayRef<const InteractionsOfType> plist,
  74                      bool                                    bFullDih)
  75 {
  76     /* This dihp is a DIRTY patch because the dih-types do not use
  77      * all four atoms to determine the type.
  78      */
  79     const int dihp[2][2] = { { 1, 2 }, { 0, 3 } };
  80     int       nral, nrfp;
  81     bool      bDih = false, bSwapParity;
  82
  83     const InteractionsOfType* bt = &(plist[ftype]);
  84
  85     if (bt->size() == 0)
  86     {
  87         return;
  88     }
  89
  90     int f = 0;
  91     switch (ftype)
  92     {
  93         case F_G96ANGLES: // Intended to fall through
  94         case F_G96BONDS: f = 1; break;
  95         case F_MORSE: f = 2; break;
  96         case F_CUBICBONDS: f = 3; break;
  97         case F_CONNBONDS: f = 4; break;
  98         case F_HARMONIC: f = 5; break;
  99         case F_CROSS_BOND_ANGLES: f = 2; break;
 100         case F_CROSS_BOND_BONDS: f = 3; break;
 101         case F_UREY_BRADLEY: f = 4; break;
 102         case F_PDIHS:  // Intended to fall through
 103         case F_RBDIHS: // Intended to fall through
 104         case F_FOURDIHS: bDih = TRUE; break;
 105         case F_IDIHS:
 106             f    = 1;
 107             bDih = TRUE;
 108             break;
 109         case F_CONSTRNC: // Intended to fall through
 110         case F_VSITE3FD: f = 1; break;
 111         case F_VSITE3FAD: f = 2; break;
 112         case F_VSITE3OUT: f = 3; break;
 113         case F_VSITE4FDN: // Intended to fall through
 114         case F_CMAP: f = 1; break;
 115
 116         default: bDih = FALSE;
 117     }
 118     if (bFullDih)
 119     {
 120         bDih = FALSE;
 121     }
 122     if (fsubtype)
 123     {
 124         f = fsubtype - 1;
 125     }
 126
 127     nral = NRAL(ftype);
 128     nrfp = NRFP(ftype);
 129
 130     /* header */
 131     fprintf(out, "[ %s ]\n", dir2str(d));
 132     fprintf(out, "; ");
 133     if (!bDih)
 134     {
 135         fprintf(out, "%3s  %4s", "ai", "aj");
 136         for (int j = 2; (j < nral); j++)
 137         {
 138             fprintf(out, "  %3c%c", 'a', 'i' + j);
 139         }
 140     }
 141     else
 142     {
 143         for (int j = 0; (j < 2); j++)
 144         {
 145             fprintf(out, "%3c%c", 'a', 'i' + dihp[f][j]);
 146         }
 147     }
 148
 149     fprintf(out, " funct");
 150     for (int j = 0; (j < nrfp); j++)
 151     {
 152         fprintf(out, " %12c%1d", 'c', j);
 153     }
 154     fprintf(out, "\n");
 155
 156     /* print bondtypes */
 157     for (const auto& parm : bt->interactionTypes)
 158     {
 159         bSwapParity                    = (parm.c0() == NOTSET) && (parm.c1() == -1);
 160         gmx::ArrayRef<const int> atoms = parm.atoms();
 161         if (!bDih)
 162         {
 163             for (int j = 0; (j < nral); j++)
 164             {
 165                 fprintf(out, "%5s ", at->atomNameFromAtomType(atoms[j]));
 166             }
 167         }
 168         else
 169         {
 170             for (int j = 0; (j < 2); j++)
 171             {
 172                 fprintf(out, "%5s ", at->atomNameFromAtomType(atoms[dihp[f][j]]));
 173             }
 174         }
 175         fprintf(out, "%5d ", bSwapParity ? -f - 1 : f + 1);
 176
 177         if (!parm.interactionTypeName().empty())
 178         {
 179             fprintf(out, "   %s", parm.interactionTypeName().c_str());
 180         }
 181         else
 182         {
 183             gmx::ArrayRef<const real> forceParam = parm.forceParam();
 184             for (int j = 0; (j < nrfp) && (forceParam[j] != NOTSET); j++)
 185             {
 186                 fprintf(out, "%13.6e ", forceParam[j]);
 187             }
 188         }
 189
 190         fprintf(out, "\n");
 191     }
 192     fprintf(out, "\n");
 193     fflush(out);
 194 }
 195
 196 void print_excl(FILE* out, int natoms, t_excls excls[])
 197 {
 198     int  i;
 199     bool have_excl;
 200     int  j;
 201
 202     have_excl = FALSE;
 203     for (i = 0; i < natoms && !have_excl; i++)
 204     {
 205         have_excl = (excls[i].nr > 0);
 206     }
 207
 208     if (have_excl)
 209     {
 210         fprintf(out, "[ %s ]\n", dir2str(Directive::d_exclusions));
 211         fprintf(out, "; %4s    %s\n", "i", "excluded from i");
 212         for (i = 0; i < natoms; i++)
 213         {
 214             if (excls[i].nr > 0)
 215             {
 216                 fprintf(out, "%6d ", i + 1);
 217                 for (j = 0; j < excls[i].nr; j++)
 218                 {
 219                     fprintf(out, " %5d", excls[i].e[j] + 1);
 220                 }
 221                 fprintf(out, "\n");
 222             }
 223         }
 224         fprintf(out, "\n");
 225         fflush(out);
 226     }
 227 }
 228
 229 static double get_residue_charge(const t_atoms* atoms, int at)
 230 {
 231     int    ri;
 232     double q;
 233
 234     ri = atoms->atom[at].resind;
 235     q  = 0;
 236     while (at < atoms->nr && atoms->atom[at].resind == ri)
 237     {
 238         q += atoms->atom[at].q;
 239         at++;
 240     }
 241
 242     return q;
 243 }
 244
 245 void print_atoms(FILE* out, PreprocessingAtomTypes* atype, t_atoms* at, int* cgnr, bool bRTPresname)
 246 {
 247     int         i, ri;
 248     int         tpA, tpB;
 249     const char* as;
 250     const char *tpnmA, *tpnmB;
 251     double      qres, qtot;
 252
 253     as = dir2str(Directive::d_atoms);
 254     fprintf(out, "[ %s ]\n", as);
 255     fprintf(out,
 256             "; %4s %10s %6s %7s%6s %6s %10s %10s %6s %10s %10s\n",
 257             "nr",
 258             "type",
 259             "resnr",
 260             "residue",
 261             "atom",
 262             "cgnr",
 263             "charge",
 264             "mass",
 265             "typeB",
 266             "chargeB",
 267             "massB");
 268
 269     qtot = 0;
 270
 271     if (at->nres)
 272     {
 273         /* if the information is present... */
 274         for (i = 0; (i < at->nr); i++)
 275         {
 276             ri = at->atom[i].resind;
 277             if ((i == 0 || ri != at->atom[i - 1].resind) && at->resinfo[ri].rtp != nullptr)
 278             {
 279                 qres = get_residue_charge(at, i);
 280                 fprintf(out,
 281                         "; residue %3d %-3s rtp %-4s q ",
 282                         at->resinfo[ri].nr,
 283                         *at->resinfo[ri].name,
 284                         *at->resinfo[ri].rtp);
 285                 if (fabs(qres) < 0.001)
 286                 {
 287                     fprintf(out, " %s", "0.0");
 288                 }
 289                 else
 290                 {
 291                     fprintf(out, "%+3.1f", qres);
 292                 }
 293                 fprintf(out, "\n");
 294             }
 295             tpA = at->atom[i].type;
 296             if ((tpnmA = atype->atomNameFromAtomType(tpA)) == nullptr)
 297             {
 298                 gmx_fatal(FARGS, "tpA = %d, i= %d in print_atoms", tpA, i);
 299             }
 300
 301             /* This is true by construction, but static analysers don't know */
 302             GMX_ASSERT(!bRTPresname || at->resinfo[at->atom[i].resind].rtp,
 303                        "-rtpres did not have residue name available");
 304             fprintf(out,
 305                     "%6d %10s %6d%c %5s %6s %6d %10g %10g",
 306                     i + 1,
 307                     tpnmA,
 308                     at->resinfo[ri].nr,
 309                     at->resinfo[ri].ic,
 310                     bRTPresname ? *(at->resinfo[at->atom[i].resind].rtp)
 311                                 : *(at->resinfo[at->atom[i].resind].name),
 312                     *(at->atomname[i]),
 313                     cgnr[i],
 314                     at->atom[i].q,
 315                     at->atom[i].m);
 316             if (PERTURBED(at->atom[i]))
 317             {
 318                 tpB = at->atom[i].typeB;
 319                 if ((tpnmB = atype->atomNameFromAtomType(tpB)) == nullptr)
 320                 {
 321                     gmx_fatal(FARGS, "tpB = %d, i= %d in print_atoms", tpB, i);
 322                 }
 323                 fprintf(out, " %6s %10g %10g", tpnmB, at->atom[i].qB, at->atom[i].mB);
 324             }
 325             // Accumulate the total charge to help troubleshoot issues.
 326             qtot += static_cast<double>(at->atom[i].q);
 327             // Round it to zero if it is close to zero, because
 328             // printing -9.34e-5 confuses users.
 329             if (fabs(qtot) < 0.0001)
 330             {
 331                 qtot = 0;
 332             }
 333             // Write the total charge for the last atom of the system
 334             // and/or residue, because generally that's where it is
 335             // expected to be an integer.
 336             if (i == at->nr - 1 || ri != at->atom[i + 1].resind)
 337             {
 338                 fprintf(out, "   ; qtot %.4g\n", qtot);
 339             }
 340             else
 341             {
 342                 fputs("\n", out);
 343             }
 344         }
 345     }
 346     fprintf(out, "\n");
 347     fflush(out);
 348 }
 349
 350 void print_bondeds(FILE* out, int natoms, Directive d, int ftype, int fsubtype, gmx::ArrayRef<const InteractionsOfType> plist)
 351 {
 352     t_symtab stab;
 353     t_atom*  a;
 354
 355     PreprocessingAtomTypes atype;
 356     snew(a, 1);
 357     open_symtab(&stab);
 358     for (int i = 0; (i < natoms); i++)
 359     {
 360         char buf[12];
 361         sprintf(buf, "%4d", (i + 1));
 362         atype.addType(&stab, *a, buf, InteractionOfType({}, {}), 0, 0);
 363     }
 364     print_bt(out, d, &atype, ftype, fsubtype, plist, TRUE);
 365
 366     done_symtab(&stab);
 367     sfree(a);
 368 }