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