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