src/gromacs/gmxpreprocess/topdirs.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) 2013,2014,2015,2017,2018 by the GROMACS development team.
   7  * Copyright (c) 2019,2020,2021, 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 "topdirs.h"
  41
  42 #include <cstdarg>
  43 #include <cstdio>
  44
  45 #include <algorithm>
  46
  47 #include "gromacs/topology/idef.h"
  48 #include "gromacs/utility/cstringutil.h"
  49 #include "gromacs/utility/enumerationhelpers.h"
  50 #include "gromacs/utility/fatalerror.h"
  51 #include "gromacs/utility/smalloc.h"
  52 #include "gromacs/utility/stringtoenumvalueconverter.h"
  53
  54 const char* enumValueToString(Directive d)
  55 {
  56     /* Must correspond to the Directive enum in topdirs.h */
  57     static constexpr gmx::EnumerationArray<Directive, const char*> directiveNames = {
  58         "defaults",
  59         "atomtypes",
  60         "bondtypes",
  61         "constrainttypes",
  62         "pairtypes",
  63         "angletypes",
  64         "dihedraltypes",
  65         "nonbond_params",
  66         "implicit_genborn_params",
  67         "implicit_surface_params",
  68         "cmaptypes",
  69         /* All the directives above can not appear after moleculetype */
  70         "moleculetype",
  71         "atoms",
  72         "virtual_sites1",
  73         "virtual_sites2",
  74         "virtual_sites3",
  75         "virtual_sites4",
  76         "virtual_sitesn",
  77         "bonds",
  78         "exclusions",
  79         "pairs",
  80         "pairs_nb",
  81         "angles",
  82         "dihedrals",
  83         "constraints",
  84         "settles",
  85         "polarization",
  86         "water_polarization",
  87         "thole_polarization",
  88         "system",
  89         "molecules",
  90         "position_restraints",
  91         "angle_restraints",
  92         "angle_restraints_z",
  93         "distance_restraints",
  94         "orientation_restraints",
  95         "dihedral_restraints",
  96         "cmap",
  97         "intermolecular_interactions",
  98         "maxdirs",
  99         "invalid",
 100         "none"
 101     };
 102     return directiveNames[d];
 103 }
 104
 105 int ifunc_index(Directive d, int type)
 106 {
 107     switch (d)
 108     {
 109         case Directive::d_bondtypes:
 110         case Directive::d_bonds:
 111             switch (type)
 112             {
 113                 case 1: return F_BONDS;
 114                 case 2: return F_G96BONDS;
 115                 case 3: return F_MORSE;
 116                 case 4: return F_CUBICBONDS;
 117                 case 5: return F_CONNBONDS;
 118                 case 6: return F_HARMONIC;
 119                 case 7: return F_FENEBONDS;
 120                 case 8: return F_TABBONDS;
 121                 case 9: return F_TABBONDSNC;
 122                 case 10: return F_RESTRBONDS;
 123                 default: gmx_fatal(FARGS, "Invalid bond type %d", type);
 124             }
 125         case Directive::d_angles:
 126         case Directive::d_angletypes:
 127             switch (type)
 128             {
 129                 case 1: return F_ANGLES;
 130                 case 2: return F_G96ANGLES;
 131                 case 3: return F_CROSS_BOND_BONDS;
 132                 case 4: return F_CROSS_BOND_ANGLES;
 133                 case 5: return F_UREY_BRADLEY;
 134                 case 6: return F_QUARTIC_ANGLES;
 135                 case 8: return F_TABANGLES;
 136                 case 9: return F_LINEAR_ANGLES;
 137                 case 10: return F_RESTRANGLES;
 138                 default: gmx_fatal(FARGS, "Invalid angle type %d", type);
 139             }
 140         case Directive::d_pairs:
 141         case Directive::d_pairtypes:
 142             if (type == 1 || (d == Directive::d_pairtypes && type == 2))
 143             {
 144                 return F_LJ14;
 145             }
 146             else if (type == 2)
 147             {
 148                 return F_LJC14_Q;
 149             }
 150             else
 151             {
 152                 gmx_fatal(FARGS, "Invalid pairs type %d", type);
 153             }
 154         case Directive::d_pairs_nb: return F_LJC_PAIRS_NB;
 155         case Directive::d_dihedrals:
 156         case Directive::d_dihedraltypes:
 157             switch (type)
 158             {
 159                 case 1: return F_PDIHS;
 160                 case 2: return F_IDIHS;
 161                 case 3: return F_RBDIHS;
 162                 case 4: return F_PIDIHS;
 163                 case 5: return F_FOURDIHS;
 164                 case 8: return F_TABDIHS;
 165                 case 9:
 166                     return F_PDIHS; /* proper dihedrals where we allow multiple terms over single bond */
 167                 case 10: return F_RESTRDIHS;
 168                 case 11: return F_CBTDIHS;
 169                 default: gmx_fatal(FARGS, "Invalid dihedral type %d", type);
 170             }
 171         case Directive::d_cmaptypes:
 172         case Directive::d_cmap: return F_CMAP;
 173
 174         case Directive::d_nonbond_params:
 175             if (type == 1)
 176             {
 177                 return F_LJ;
 178             }
 179             else
 180             {
 181                 return F_BHAM;
 182             }
 183         case Directive::d_vsites1:
 184             if (type == 1)
 185             {
 186                 return F_VSITE1;
 187             }
 188             else
 189             {
 190                 gmx_fatal(FARGS, "Invalid vsites1 type %d", type);
 191             }
 192         case Directive::d_vsites2:
 193             switch (type)
 194             {
 195                 case 1: return F_VSITE2;
 196                 case 2: return F_VSITE2FD;
 197                 default: gmx_fatal(FARGS, "Invalid vsites2 type %d", type);
 198             }
 199         case Directive::d_vsites3:
 200             switch (type)
 201             {
 202                 case 1: return F_VSITE3;
 203                 case 2: return F_VSITE3FD;
 204                 case 3: return F_VSITE3FAD;
 205                 case 4: return F_VSITE3OUT;
 206                 default: gmx_fatal(FARGS, "Invalid vsites3 type %d", type);
 207             }
 208         case Directive::d_vsites4:
 209             switch (type)
 210             {
 211                 case 1: return F_VSITE4FD;
 212                 case 2: return F_VSITE4FDN;
 213                 default: gmx_fatal(FARGS, "Invalid vsites4 type %d", type);
 214             }
 215         case Directive::d_vsitesn: return F_VSITEN;
 216         case Directive::d_constraints:
 217         case Directive::d_constrainttypes:
 218             switch (type)
 219             {
 220                 case 1: return F_CONSTR;
 221                 case 2: return F_CONSTRNC;
 222                 default: gmx_fatal(FARGS, "Invalid constraints type %d", type);
 223             }
 224         case Directive::d_settles: return F_SETTLE;
 225         case Directive::d_position_restraints:
 226             switch (type)
 227             {
 228                 case 1: return F_POSRES;
 229                 case 2: return F_FBPOSRES;
 230                 default: gmx_fatal(FARGS, "Invalid position restraint type %d", type);
 231             }
 232         case Directive::d_polarization:
 233             switch (type)
 234             {
 235                 case 1: return F_POLARIZATION;
 236                 case 2: return F_ANHARM_POL;
 237                 default: gmx_fatal(FARGS, "Invalid polarization type %d", type);
 238             }
 239         case Directive::d_thole_polarization: return F_THOLE_POL;
 240         case Directive::d_water_polarization: return F_WATER_POL;
 241         case Directive::d_angle_restraints: return F_ANGRES;
 242         case Directive::d_angle_restraints_z: return F_ANGRESZ;
 243         case Directive::d_distance_restraints: return F_DISRES;
 244         case Directive::d_orientation_restraints: return F_ORIRES;
 245         case Directive::d_dihedral_restraints: return F_DIHRES;
 246         default:
 247             gmx_fatal(FARGS, "invalid directive %s in ifunc_index (%s:%d)", enumValueToString(d), __FILE__, __LINE__);
 248     }
 249 }
 250
 251 enum class DeprecatedDirectives : int
 252 {
 253     d_dummies1,
 254     d_dummies2,
 255     d_dummies3,
 256     d_dummies4,
 257     d_dummiesn,
 258     Count
 259 };
 260
 261 static const char* enumValueToString(DeprecatedDirectives d)
 262 {
 263     static constexpr gmx::EnumerationArray<DeprecatedDirectives, const char*> directiveNames = {
 264         "dummies1", "dummies2", "dummies3", "dummies4", "dummiesn"
 265     };
 266     return directiveNames[d];
 267 }
 268
 269 Directive str2dir(char* dstr)
 270 {
 271     static const gmx::StringToEnumValueConverter<Directive, enumValueToString, gmx::StringCompareType::CaseAndDashInsensitive> s_converter;
 272
 273     if (std::optional<Directive> d = s_converter.valueFrom(dstr); d.has_value())
 274     {
 275         return d.value();
 276     }
 277     // Also handle deprecated directives that have modern replacements, like
 278     // "dummies*" -> "virtual_sites*"
 279
 280     static const gmx::StringToEnumValueConverter<DeprecatedDirectives, enumValueToString, gmx::StringCompareType::CaseAndDashInsensitive>
 281             s_converterForDeprecated;
 282
 283     if (std::optional<DeprecatedDirectives> d = s_converterForDeprecated.valueFrom(dstr); d.has_value())
 284     {
 285         static constexpr gmx::EnumerationArray<DeprecatedDirectives, Directive> s_deprecatedDirectiveToDirective = {
 286             Directive::d_vsites1, Directive::d_vsites2, Directive::d_vsites3,
 287             Directive::d_vsites4, Directive::d_vsitesn,
 288         };
 289         return s_deprecatedDirectiveToDirective[d.value()];
 290     }
 291     return Directive::d_invalid;
 292 }
 293
 294 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
 295 static gmx::EnumerationArray<Directive, Directive*> necessary = { { nullptr } };
 296
 297 static void set_nec(Directive** n, ...)
 298 /* Must always have at least one extra argument */
 299 {
 300     va_list   ap;
 301     int       ind = 0;
 302     Directive d;
 303
 304     va_start(ap, n);
 305     do
 306     {
 307         d = static_cast<Directive>(va_arg(ap, int));
 308         srenew(*n, ++ind);
 309         (*n)[ind - 1] = d;
 310     } while (d != Directive::d_none);
 311     va_end(ap);
 312 }
 313
 314 void DS_Init(DirStack** DS)
 315 {
 316     if (necessary[0] == nullptr)
 317     {
 318         set_nec(&(necessary[Directive::d_defaults]), Directive::d_none);
 319         set_nec(&(necessary[Directive::d_atomtypes]), Directive::d_defaults, Directive::d_none);
 320         set_nec(&(necessary[Directive::d_bondtypes]), Directive::d_atomtypes, Directive::d_none);
 321         set_nec(&(necessary[Directive::d_constrainttypes]), Directive::d_atomtypes, Directive::d_none);
 322         set_nec(&(necessary[Directive::d_pairtypes]), Directive::d_atomtypes, Directive::d_none);
 323         set_nec(&(necessary[Directive::d_angletypes]), Directive::d_atomtypes, Directive::d_none);
 324         set_nec(&(necessary[Directive::d_dihedraltypes]), Directive::d_atomtypes, Directive::d_none);
 325         set_nec(&(necessary[Directive::d_nonbond_params]), Directive::d_atomtypes, Directive::d_none);
 326         // Note that the content of the next two directives are
 327         // ignored, but if grompp reads them in old force field files,
 328         // it still needs to understand that they are in a valid place
 329         // in the .top structure. It doesn't have to require them to
 330         // be in the same place that was valid in old versions (ie. child
 331         // directive of [atomtypes]) but any relevant case will
 332         // satisfy that.
 333         set_nec(&(necessary[Directive::d_implicit_genborn_params]), Directive::d_atomtypes, Directive::d_none);
 334         set_nec(&(necessary[Directive::d_implicit_surface_params]), Directive::d_atomtypes, Directive::d_none);
 335         set_nec(&(necessary[Directive::d_cmaptypes]), Directive::d_atomtypes, Directive::d_none);
 336         set_nec(&(necessary[Directive::d_moleculetype]), Directive::d_atomtypes, Directive::d_none);
 337         set_nec(&(necessary[Directive::d_atoms]), Directive::d_moleculetype, Directive::d_none);
 338         set_nec(&(necessary[Directive::d_vsites1]), Directive::d_atoms, Directive::d_none);
 339         set_nec(&(necessary[Directive::d_vsites2]), Directive::d_atoms, Directive::d_none);
 340         set_nec(&(necessary[Directive::d_vsites3]), Directive::d_atoms, Directive::d_none);
 341         set_nec(&(necessary[Directive::d_vsites4]), Directive::d_atoms, Directive::d_none);
 342         set_nec(&(necessary[Directive::d_vsitesn]), Directive::d_atoms, Directive::d_none);
 343         set_nec(&(necessary[Directive::d_bonds]), Directive::d_atoms, Directive::d_none);
 344         set_nec(&(necessary[Directive::d_exclusions]),
 345                 Directive::d_bonds,
 346                 Directive::d_constraints,
 347                 Directive::d_settles,
 348                 Directive::d_none);
 349         set_nec(&(necessary[Directive::d_pairs]), Directive::d_atoms, Directive::d_none);
 350         set_nec(&(necessary[Directive::d_pairs_nb]), Directive::d_atoms, Directive::d_none);
 351         set_nec(&(necessary[Directive::d_angles]), Directive::d_atoms, Directive::d_none);
 352         set_nec(&(necessary[Directive::d_polarization]), Directive::d_atoms, Directive::d_none);
 353         set_nec(&(necessary[Directive::d_water_polarization]), Directive::d_atoms, Directive::d_none);
 354         set_nec(&(necessary[Directive::d_thole_polarization]), Directive::d_atoms, Directive::d_none);
 355         set_nec(&(necessary[Directive::d_dihedrals]), Directive::d_atoms, Directive::d_none);
 356         set_nec(&(necessary[Directive::d_constraints]), Directive::d_atoms, Directive::d_none);
 357         set_nec(&(necessary[Directive::d_settles]), Directive::d_atoms, Directive::d_none);
 358         set_nec(&(necessary[Directive::d_system]), Directive::d_moleculetype, Directive::d_none);
 359         set_nec(&(necessary[Directive::d_molecules]), Directive::d_system, Directive::d_none);
 360         set_nec(&(necessary[Directive::d_position_restraints]), Directive::d_atoms, Directive::d_none);
 361         set_nec(&(necessary[Directive::d_angle_restraints]), Directive::d_atoms, Directive::d_none);
 362         set_nec(&(necessary[Directive::d_angle_restraints_z]), Directive::d_atoms, Directive::d_none);
 363         set_nec(&(necessary[Directive::d_distance_restraints]), Directive::d_atoms, Directive::d_none);
 364         set_nec(&(necessary[Directive::d_orientation_restraints]), Directive::d_atoms, Directive::d_none);
 365         set_nec(&(necessary[Directive::d_dihedral_restraints]), Directive::d_atoms, Directive::d_none);
 366         set_nec(&(necessary[Directive::d_cmap]), Directive::d_atoms, Directive::d_none);
 367         set_nec(&(necessary[Directive::d_intermolecular_interactions]),
 368                 Directive::d_molecules,
 369                 Directive::d_none);
 370     }
 371     *DS = nullptr;
 372 }
 373
 374 void DS_Done(DirStack** DS)
 375 {
 376     DirStack* D;
 377
 378     while (*DS != nullptr)
 379     {
 380         D   = *DS;
 381         *DS = (*DS)->prev;
 382         sfree(D);
 383     }
 384 }
 385
 386 void DS_Push(DirStack** DS, Directive d)
 387 {
 388     DirStack* D;
 389
 390     snew(D, 1);
 391     D->d    = d;
 392     D->prev = *DS;
 393     *DS     = D;
 394 }
 395
 396 int DS_Search(DirStack* DS, Directive d)
 397 {
 398     DirStack* D;
 399
 400     D = DS;
 401     while ((D != nullptr) && (D->d != d))
 402     {
 403         D = D->prev;
 404     }
 405
 406     return static_cast<int>(D != nullptr);
 407 }
 408
 409 int DS_Check_Order(DirStack* DS, Directive d)
 410 {
 411     Directive d0;
 412     int       i = 0;
 413
 414     /* Check if parameter definitions appear after a moleculetype directive */
 415     if (d < Directive::d_moleculetype && DS_Search(DS, Directive::d_moleculetype))
 416     {
 417         return FALSE;
 418     }
 419
 420     /* Check if all the necessary directives have appeared before directive d */
 421     if (necessary[d][0] == Directive::d_none)
 422     {
 423         return TRUE;
 424     }
 425     else
 426     {
 427         do
 428         {
 429             d0 = necessary[d][i++];
 430             if (DS_Search(DS, d0))
 431             {
 432                 return TRUE;
 433             }
 434         } while (d0 != Directive::d_none);
 435     }
 436     return FALSE;
 437 }