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