src/gromacs/gmxpreprocess/topdirs.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <stdio.h>
  40 #include <stdarg.h>
  41
  42 #include "sysstuff.h"
  43 #include "smalloc.h"
  44 #include "macros.h"
  45 #include "string2.h"
  46 #include "gmx_fatal.h"
  47 #include "topdirs.h"
  48
  49 /* Must correspond to the directive enum in grompp.h */
  50 static const char *directive_names[d_maxdir+1] = {
  51     "defaults",
  52     "atomtypes",
  53     "bondtypes",
  54     "constrainttypes",
  55     "pairtypes",
  56     "angletypes",
  57     "dihedraltypes",
  58     "nonbond_params",
  59     "implicit_genborn_params",
  60     "implicit_surface_params",
  61     "cmaptypes",
  62     /* All the directives above can not appear after moleculetype */
  63     "moleculetype",
  64     "atoms",
  65     "virtual_sites2",
  66     "virtual_sites3",
  67     "virtual_sites4",
  68     "virtual_sitesn",
  69     "bonds",
  70     "exclusions",
  71     "pairs",
  72     "pairs_nb",
  73     "angles",
  74     "dihedrals",
  75     "constraints",
  76     "settles",
  77     "polarization",
  78     "water_polarization",
  79     "thole_polarization",
  80     "system",
  81     "molecules",
  82     "position_restraints",
  83     "angle_restraints",
  84     "angle_restraints_z",
  85     "distance_restraints",
  86     "orientation_restraints",
  87     "dihedral_restraints",
  88     "cmap",
  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                     break;
 123             }
 124             break;
 125         case d_angles:
 126         case d_angletypes:
 127             switch (type)
 128             {
 129                 case 1:
 130                     return F_ANGLES;
 131                 case 2:
 132                     return F_G96ANGLES;
 133                 case 3:
 134                     return F_CROSS_BOND_BONDS;
 135                 case 4:
 136                     return F_CROSS_BOND_ANGLES;
 137                 case 5:
 138                     return F_UREY_BRADLEY;
 139                 case 6:
 140                     return F_QUARTIC_ANGLES;
 141                 case 8:
 142                     return F_TABANGLES;
 143                 case 9:
 144                     return F_LINEAR_ANGLES;
 145                 default:
 146                     gmx_fatal(FARGS, "Invalid angle type %d", type);
 147                     break;
 148             }
 149             break;
 150         case d_pairs:
 151         case d_pairtypes:
 152             if (type == 1 || (d == d_pairtypes && type == 2))
 153             {
 154                 return F_LJ14;
 155             }
 156             else if (type == 2)
 157             {
 158                 return F_LJC14_Q;
 159             }
 160             else
 161             {
 162                 gmx_fatal(FARGS, "Invalid pairs type %d", type);
 163             }
 164             break;
 165         case d_pairs_nb:
 166             return F_LJC_PAIRS_NB;
 167         case d_dihedrals:
 168         case d_dihedraltypes:
 169             switch (type)
 170             {
 171                 case 1:
 172                     return F_PDIHS;
 173                 case 2:
 174                     return F_IDIHS;
 175                 case 3:
 176                     return F_RBDIHS;
 177                 case 4:
 178                     return F_PIDIHS;
 179                 case 5:
 180                     return F_FOURDIHS;
 181                 case 8:
 182                     return F_TABDIHS;
 183                 case 9:
 184                     return F_PDIHS; /* proper dihedrals where we allow multiple terms over single bond */
 185                 default:
 186                     gmx_fatal(FARGS, "Invalid dihedral type %d", type);
 187             }
 188             break;
 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             break;
 219         case d_vsites4:
 220             switch (type)
 221             {
 222                 case 1:
 223                     return F_VSITE4FD;
 224                 case 2:
 225                     return F_VSITE4FDN;
 226                 default:
 227                     gmx_fatal(FARGS, "Invalid vsites4 type %d", type);
 228             }
 229             break;
 230         case d_vsitesn:
 231             return F_VSITEN;
 232         case d_constraints:
 233         case d_constrainttypes:
 234             switch (type)
 235             {
 236                 case 1:
 237                     return F_CONSTR;
 238                 case 2:
 239                     return F_CONSTRNC;
 240                 default:
 241                     gmx_fatal(FARGS, "Invalid constraints type %d", type);
 242             }
 243             break;
 244         case d_settles:
 245             return F_SETTLE;
 246         case d_position_restraints:
 247             switch (type)
 248             {
 249                 case 1:
 250                     return F_POSRES;
 251                 case 2:
 252                     return F_FBPOSRES;
 253                     break;
 254                 default:
 255                     gmx_fatal(FARGS, "Invalid position restraint type %d", type);
 256             }
 257             break;
 258         case d_polarization:
 259             switch (type)
 260             {
 261                 case 1:
 262                     return F_POLARIZATION;
 263                 case 2:
 264                     return F_ANHARM_POL;
 265                 default:
 266                     gmx_fatal(FARGS, "Invalid polarization type %d", type);
 267             }
 268             break;
 269         case d_thole_polarization:
 270             return F_THOLE_POL;
 271         case d_water_polarization:
 272             return F_WATER_POL;
 273         case d_angle_restraints:
 274             return F_ANGRES;
 275         case d_angle_restraints_z:
 276             return F_ANGRESZ;
 277         case d_distance_restraints:
 278             return F_DISRES;
 279         case d_orientation_restraints:
 280             return F_ORIRES;
 281         case d_dihedral_restraints:
 282             return F_DIHRES;
 283         default:
 284             gmx_fatal(FARGS, "invalid directive %s in ifunc_index (%s:%s)",
 285                       dir2str(d), __FILE__, __LINE__);
 286     }
 287     return -1;
 288 }
 289
 290 const char *dir2str (directive d)
 291 {
 292     if (d < d_maxdir)
 293     {
 294         return directive_names[d];
 295     }
 296     else
 297     {
 298         return directive_names[d_maxdir];
 299     }
 300 }
 301
 302 directive str2dir (char *dstr)
 303 {
 304     int  d;
 305     char buf[STRLEN], *ptr;
 306
 307     /* Hack to be able to read old topologies */
 308     if (gmx_strncasecmp_min(dstr, "dummies", 7) == 0)
 309     {
 310         sprintf(buf, "virtual_sites%s", dstr+7);
 311         ptr = buf;
 312     }
 313     else
 314     {
 315         ptr = dstr;
 316     }
 317
 318     for (d = 0; (d < d_maxdir); d++)
 319     {
 320         if (gmx_strcasecmp_min(ptr, dir2str((directive)d)) == 0)
 321         {
 322             return (directive)d;
 323         }
 324     }
 325
 326     return d_invalid;
 327 }
 328
 329 static directive **necessary = NULL;
 330
 331 static void set_nec(directive **n, ...)
 332 /* Must always have at least one extra argument */
 333 {
 334     va_list   ap;
 335     int       ind = 0;
 336     directive d;
 337
 338     va_start(ap, n);
 339     do
 340     {
 341         d = (directive)va_arg(ap, int);
 342         srenew(*n, ++ind);
 343         (*n)[ind-1] = d;
 344     }
 345     while (d != d_none);
 346     va_end(ap);
 347 }
 348
 349 void DS_Init(DirStack **DS)
 350 {
 351     if (necessary == NULL)
 352     {
 353         int i;
 354
 355         snew(necessary, d_maxdir);
 356         set_nec(&(necessary[d_defaults]), d_none);
 357         set_nec(&(necessary[d_atomtypes]), d_defaults, d_none);
 358         set_nec(&(necessary[d_bondtypes]), d_atomtypes, d_none);
 359         set_nec(&(necessary[d_constrainttypes]), d_atomtypes, d_none);
 360         set_nec(&(necessary[d_pairtypes]), d_atomtypes, d_none);
 361         set_nec(&(necessary[d_angletypes]), d_atomtypes, d_none);
 362         set_nec(&(necessary[d_dihedraltypes]), d_atomtypes, d_none);
 363         set_nec(&(necessary[d_nonbond_params]), d_atomtypes, d_none);
 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
 394         for (i = 0; (i < d_maxdir); i++)
 395         {
 396             if (debug)
 397             {
 398                 fprintf(debug, "%20s:  ", dir2str((directive)i));
 399             }
 400             if (necessary[i])
 401             {
 402                 directive d;
 403                 int       j = 0;
 404                 do
 405                 {
 406                     d = necessary[i][j++];
 407                     if (debug)
 408                     {
 409                         fprintf(debug, "%20s  ", dir2str(d));
 410                     }
 411                 }
 412                 while (d != d_none);
 413             }
 414             if (debug)
 415             {
 416                 fprintf(debug, "\n");
 417             }
 418         }
 419     }
 420     *DS = NULL;
 421
 422 }
 423
 424 void DS_Done (DirStack **DS)
 425 {
 426     DirStack *D;
 427
 428     while (*DS != NULL)
 429     {
 430         D   = *DS;
 431         *DS = (*DS)->prev;
 432         sfree (D);
 433     }
 434 }
 435
 436 void DS_Push (DirStack **DS, directive d)
 437 {
 438     DirStack *D;
 439
 440     snew(D, 1);
 441     D->d    = d;
 442     D->prev = *DS;
 443     *DS     = D;
 444 }
 445
 446 int DS_Search(DirStack *DS, directive d)
 447 {
 448     DirStack *D;
 449
 450     D = DS;
 451     while ((D != NULL) && (D->d != d))
 452     {
 453         D = D->prev;
 454     }
 455
 456     return (D != NULL);
 457 }
 458
 459 int DS_Check_Order(DirStack *DS, directive d)
 460 {
 461     directive d0;
 462     int       i = 0;
 463
 464     /* Check if parameter definitions appear after a moleculetype directive */
 465     if (d < d_moleculetype && DS_Search(DS, d_moleculetype))
 466     {
 467         return FALSE;
 468     }
 469
 470     /* Check if all the necessary directives have appeared before directive d */
 471     if (necessary[d][0] == d_none)
 472     {
 473         return TRUE;
 474     }
 475     else
 476     {
 477         do
 478         {
 479             d0 = necessary[d][i++];
 480             if (DS_Search(DS, d0))
 481             {
 482                 return TRUE;
 483             }
 484         }
 485         while (d0 != d_none);
 486     }
 487     return FALSE;
 488 }