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