src/gromacs/gmxpreprocess/topdirs.c

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