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