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