src/kernel/topdirs.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <stdio.h>
  43 #include <stdarg.h>
  44
  45 #include "sysstuff.h"
  46 #include "smalloc.h"
  47 #include "macros.h"
  48 #include "string2.h"
  49 #include "gmx_fatal.h"
  50 #include "topdirs.h"
  51
  52 int ifunc_index(directive d,int type)
  53 {
  54   switch (d) {
  55   case d_bondtypes:
  56   case d_bonds:
  57     switch(type) {
  58     case 1:
  59       return F_BONDS;
  60     case 2:
  61       return F_G96BONDS;
  62     case 3:
  63       return F_MORSE;
  64     case 4:
  65       return F_CUBICBONDS;
  66     case 5:
  67       return F_CONNBONDS;
  68     case 6:
  69       return F_HARMONIC;
  70     case 7:
  71       return F_FENEBONDS;
  72     case 8:
  73       return F_TABBONDS;
  74     case 9:
  75       return F_TABBONDSNC;
  76     case 10:
  77       return F_RESTRBONDS;
  78     default:
  79       gmx_fatal(FARGS,"Invalid bond type %d",type);
  80       break;
  81     }
  82     break;
  83   case d_angles:
  84   case d_angletypes:
  85     switch (type) {
  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       return F_LJ14;
 111     else if (type == 2)
 112       return F_LJC14_Q;
 113     else
 114       gmx_fatal(FARGS,"Invalid pairs type %d",type);
 115     break;
 116   case d_pairs_nb:
 117     return F_LJC_PAIRS_NB;
 118   case d_dihedrals:
 119   case d_dihedraltypes:
 120     switch (type) {
 121     case 1:
 122       return F_PDIHS;
 123     case 2:
 124       return F_IDIHS;
 125     case 3:
 126       return F_RBDIHS;
 127     case 4:
 128       return F_PIDIHS;
 129     case 5:
 130       return F_FOURDIHS;
 131     case 8:
 132       return F_TABDIHS;
 133     case 9:
 134       return F_PDIHS;  /* proper dihedrals where we allow multiple terms over single bond */
 135     default:
 136       gmx_fatal(FARGS,"Invalid dihedral type %d",type);
 137     }
 138     break;
 139   case d_cmaptypes:
 140   case d_cmap:
 141     return F_CMAP;
 142
 143   case d_nonbond_params:
 144     if (type == 1)
 145       return F_LJ;
 146     else
 147       return F_BHAM;
 148   case d_vsites2:
 149     return F_VSITE2;
 150   case d_vsites3:
 151     switch (type) {
 152     case 1:
 153       return F_VSITE3;
 154     case 2:
 155       return F_VSITE3FD;
 156     case 3:
 157       return F_VSITE3FAD;
 158     case 4:
 159       return F_VSITE3OUT;
 160     default:
 161       gmx_fatal(FARGS,"Invalid vsites3 type %d",type);
 162     }
 163     break;
 164   case d_vsites4:
 165       switch (type) {
 166           case 1:
 167               return F_VSITE4FD;
 168           case 2:
 169               return F_VSITE4FDN;
 170           default:
 171               gmx_fatal(FARGS,"Invalid vsites4 type %d",type);
 172       }
 173       break;
 174   case d_vsitesn:
 175     return F_VSITEN;
 176   case d_constraints:
 177   case d_constrainttypes:
 178     switch (type) {
 179     case 1:
 180       return F_CONSTR;
 181     case 2:
 182       return F_CONSTRNC;
 183     default:
 184       gmx_fatal(FARGS,"Invalid constraints type %d",type);
 185     }
 186     break;
 187   case d_settles:
 188     return F_SETTLE;
 189   case d_position_restraints:
 190     switch (type) {
 191     case 1:
 192       return F_POSRES;
 193     case 2:
 194       gmx_fatal(FARGS,"Water polarization should now be listed under [ water_polarization ]\n");
 195       break;
 196     default:
 197       gmx_fatal(FARGS,"Invalid position restraint type %d",type);
 198     }
 199     break;
 200   case d_polarization:
 201     switch (type) {
 202     case 1:
 203       return F_POLARIZATION;
 204     case 2:
 205       return F_ANHARM_POL;
 206     default:
 207       gmx_fatal(FARGS,"Invalid polarization type %d",type);
 208     }
 209     break;
 210   case d_thole_polarization:
 211     return F_THOLE_POL;
 212   case d_water_polarization:
 213     return F_WATER_POL;
 214   case d_angle_restraints:
 215     return F_ANGRES;
 216   case d_angle_restraints_z:
 217     return F_ANGRESZ;
 218   case d_distance_restraints:
 219     return F_DISRES;
 220   case d_orientation_restraints:
 221     return F_ORIRES;
 222   case d_dihedral_restraints:
 223     return F_DIHRES;
 224   default:
 225     gmx_fatal(FARGS,"invalid directive %s in ifunc_index (%s:%s)",
 226                 dir2str(d),__FILE__,__LINE__);
 227   }
 228   return -1;
 229 }
 230
 231 const char *dir2str (directive d)
 232 {
 233   if (d < d_maxdir)
 234     return ds[d];
 235   else
 236     return ds[d_maxdir];
 237 }
 238
 239 directive str2dir (char *dstr)
 240 {
 241   int d;
 242   char buf[STRLEN],*ptr;
 243
 244   /* Hack to be able to read old topologies */
 245   if (gmx_strncasecmp_min(dstr,"dummies",7) == 0) {
 246     sprintf(buf,"virtual_sites%s",dstr+7);
 247     ptr = buf;
 248   } else {
 249     ptr = dstr;
 250   }
 251
 252   for (d=0; (d<d_maxdir); d++)
 253     if (gmx_strcasecmp_min(ptr,dir2str((directive)d)) == 0)
 254       return (directive)d;
 255
 256   return d_invalid;
 257 }
 258
 259 static directive **necessary = NULL;
 260
 261 static void set_nec(directive **n, ...)
 262 /* Must always have at least one extra argument */
 263 {
 264   va_list   ap;
 265   int       ind=0;
 266   directive d;
 267
 268   va_start(ap,n);
 269   do {
 270     d = (directive)va_arg(ap,int);
 271     srenew(*n,++ind);
 272     (*n)[ind-1]=d;
 273   } while (d != d_none);
 274   va_end(ap);
 275 }
 276
 277 void DS_Init(DirStack **DS)
 278 {
 279   if (necessary==NULL) {
 280     int i;
 281
 282     snew(necessary,d_maxdir);
 283     set_nec(&(necessary[d_defaults]),d_none);
 284     set_nec(&(necessary[d_atomtypes]),d_defaults,d_none);
 285     set_nec(&(necessary[d_bondtypes]),d_atomtypes,d_none);
 286     set_nec(&(necessary[d_constrainttypes]),d_atomtypes,d_none);
 287     set_nec(&(necessary[d_pairtypes]),d_atomtypes,d_none);
 288     set_nec(&(necessary[d_angletypes]),d_atomtypes,d_none);
 289     set_nec(&(necessary[d_dihedraltypes]),d_atomtypes,d_none);
 290     set_nec(&(necessary[d_nonbond_params]),d_atomtypes,d_none);
 291     set_nec(&(necessary[d_implicit_genborn_params]),d_atomtypes,d_none);
 292     set_nec(&(necessary[d_implicit_surface_params]),d_atomtypes,d_none);
 293     set_nec(&(necessary[d_cmaptypes]),d_atomtypes,d_none);
 294     set_nec(&(necessary[d_moleculetype]),d_atomtypes,d_none);
 295     set_nec(&(necessary[d_atoms]),d_moleculetype,d_none);
 296     set_nec(&(necessary[d_vsites2]),d_atoms,d_none);
 297     set_nec(&(necessary[d_vsites3]),d_atoms,d_none);
 298     set_nec(&(necessary[d_vsites4]),d_atoms,d_none);
 299     set_nec(&(necessary[d_vsitesn]),d_atoms,d_none);
 300     set_nec(&(necessary[d_bonds]),d_atoms,d_none);
 301     set_nec(&(necessary[d_exclusions]),d_bonds,d_constraints,d_settles,d_none);
 302     set_nec(&(necessary[d_pairs]),d_atoms,d_none);
 303     set_nec(&(necessary[d_pairs_nb]),d_atoms,d_none);
 304     set_nec(&(necessary[d_angles]),d_atoms,d_none);
 305     set_nec(&(necessary[d_polarization]),d_atoms,d_none);
 306     set_nec(&(necessary[d_water_polarization]),d_atoms,d_none);
 307     set_nec(&(necessary[d_thole_polarization]),d_atoms,d_none);
 308     set_nec(&(necessary[d_dihedrals]),d_atoms,d_none);
 309     set_nec(&(necessary[d_constraints]),d_atoms,d_none);
 310     set_nec(&(necessary[d_settles]),d_atoms,d_none);
 311     set_nec(&(necessary[d_system]),d_moleculetype,d_none);
 312     set_nec(&(necessary[d_molecules]),d_system,d_none);
 313     set_nec(&(necessary[d_position_restraints]),d_atoms,d_none);
 314     set_nec(&(necessary[d_angle_restraints]),d_atoms,d_none);
 315     set_nec(&(necessary[d_angle_restraints_z]),d_atoms,d_none);
 316     set_nec(&(necessary[d_distance_restraints]),d_atoms,d_none);
 317     set_nec(&(necessary[d_orientation_restraints]),d_atoms,d_none);
 318     set_nec(&(necessary[d_dihedral_restraints]),d_atoms,d_none);
 319     set_nec(&(necessary[d_cmap]), d_atoms, d_none);
 320
 321     for(i=0; (i<d_maxdir); i++) {
 322       if (debug)
 323         fprintf(debug,"%20s:  ",dir2str((directive)i));
 324       if(necessary[i]) {
 325         directive d;
 326         int       j=0;
 327         do {
 328           d=necessary[i][j++];
 329           if (debug)
 330             fprintf(debug,"%20s  ",dir2str(d));
 331         } while (d!=d_none);
 332       }
 333       if (debug)
 334         fprintf(debug,"\n");
 335     }
 336   }
 337   *DS = NULL;
 338
 339 }
 340
 341 void DS_Done (DirStack **DS)
 342 {
 343   DirStack *D;
 344
 345   while (*DS != NULL) {
 346     D = *DS;
 347     *DS = (*DS)->prev;
 348     sfree (D);
 349   }
 350 }
 351
 352 void DS_Push (DirStack **DS, directive d)
 353 {
 354   DirStack *D;
 355
 356   snew(D,1);
 357   D->d = d;
 358   D->prev = *DS;
 359   *DS = D;
 360 }
 361
 362 int DS_Search(DirStack *DS, directive d)
 363 {
 364   DirStack *D;
 365
 366   D = DS;
 367   while ((D != NULL) && (D->d != d))
 368     D = D->prev;
 369
 370   return (D != NULL);
 371 }
 372
 373 int DS_Check_Order(DirStack *DS,directive d)
 374 {
 375   directive d0;
 376   int       i=0;
 377
 378   /* Check if parameter definitions appear after a moleculetype directive */
 379   if (d<d_moleculetype && DS_Search(DS,d_moleculetype))
 380     return FALSE;
 381
 382   /* Check if all the necessary directives have appeared before directive d */
 383   if (necessary[d][0] == d_none)
 384     return TRUE;
 385   else {
 386     do {
 387       d0=necessary[d][i++];
 388       if (DS_Search(DS,d0))
 389         return TRUE;
 390     } while(d0!=d_none);
 391   }
 392   return FALSE;
 393 }
 394
 395
 396