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