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