GromacsTopFile.C

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "common.h"
#include "ResizeArray.h"
#include "GromacsTopFile.h"

#define JOULES_PER_CALORIE 4.184

/* A GromacsTopFile represents the information stored in a GROMACS
   topolgy file.  This is an immutable type. */

/* XXX warning: this code contains a few algorithms which run in a
   time of O(n^3) or so in the size of the topology file, since I
   haven't bothered to do any sorting.  I don't think it matters much,
   since it still manages to run on the biggest files in less than a
   second, and nothing (memory or time) depends on the total number of
   atoms in the simulation - all that matters is the number that have
   to be individually defined. */

/* GromacsTopFile initializes this to represent the data stored in the
   file <filename>, or exits on error. */
#define LINESIZE 200

/* modes */
#define UNKNOWN 0
#define ATOMS 1
#define MOLECULETYPE 2
#define MOLECULES 3
#define SYSTEM 4
#define BONDS 5
#define BONDTYPES 6
#define ANGLES 7
#define ATOMTYPES 8
#define ANGLETYPES 9
#define DIHEDRALS 10
#define DIHEDRALTYPES 11
#define DEFAULTS 12
#define NONBOND 13
#define PAIR

GromacsTopFile::GromacsTopFile(char *filename) {
  /* open the file */
  FILE *f = fopen(filename,"r");
  char buf[LINESIZE];
  char modename[20];
  int mode;
  if(f==NULL) {
    sprintf(buf,"Error opening file '%s'",filename);
    NAMD_die(buf);
  }

  /* really bad parser XXX probably just works on the files we
     happen to have REWRITE THIS SOON.  It should allow for \- line
     continuations, check for errors in the file, etc. */
  while(fgets(buf,LINESIZE-1,f)) {
    char testchar;
    int i,j;

    /* defaults */
    int nbfunc, combrule;
    char genpairs[20];
    
    /* atom buffers */
    int num, resnum, chargegp, typenum;
    char type[NAMESIZE+1], restype[NAMESIZE+1], atomname[NAMESIZE+1];
    char particletype[NAMESIZE+1];
    float charge, mass, c6, c12, junkf;

    /* moltype buffers */
    int nrexcl;
    char molname[LONGNAMESIZE+1];

    /* molInst buffers */
    int copies;
    MolInst *moleculeinstance;

    /* general atomset buffers */
    int atomi, atomj, atomk, atoml;
    char typea[NAMESIZE+1],typeb[NAMESIZE+1],
      typec[NAMESIZE+1],typed[NAMESIZE+1];
    const char *tmptypea,*tmptypeb,*tmptypec,*tmptyped;
    int funct, index;
    float c0,c1;
    
    /* bonds */
    float b0,kB,th0,kth;

    /* dihedrals */
    float c[6];
    int mult=0;

    /* check for comments */
    if(sscanf(buf," %c",&testchar)==1) {
      if(testchar == ';') continue;
    }
    else { /* this is a blank line */
      continue;
    }

    /* check for a new mode */
    if(sscanf(buf," [ %19[^] ] ]",modename)==1) {
      /* switch the mode */
      if(0==strcmp(modename,"atoms"))              mode = ATOMS;
      else if(0==strcmp(modename,"atomtypes"))     mode = ATOMTYPES;
      else if(0==strcmp(modename,"moleculetype"))  mode = MOLECULETYPE;
      else if(0==strcmp(modename,"molecules"))     mode = MOLECULES;
      else if(0==strcmp(modename,"system"))        mode = SYSTEM;
      else if(0==strcmp(modename,"bonds"))         mode = BONDS;
      else if(0==strcmp(modename,"bondtypes"))     mode = BONDTYPES;
      else if(0==strcmp(modename,"angles"))        mode = ANGLES;
      else if(0==strcmp(modename,"angletypes"))    mode = ANGLETYPES;
      else if(0==strcmp(modename,"dihedrals"))     mode = DIHEDRALS;
      else if(0==strcmp(modename,"dihedraltypes")) mode = DIHEDRALTYPES;
      else if(0==strcmp(modename,"defaults"))      mode = DEFAULTS;
      else if(0==strcmp(modename,"nonbond_params")) mode = NONBOND;
      else {    
        fprintf(stderr,"Warning: unknown mode %s\n",modename);
        mode = UNKNOWN;
      }

      continue;
    }

    /* now do the appropriate thing based on the current mode */
    switch(mode) {
    case SYSTEM:
      systemName = strdup(buf);
      break;

    case DEFAULTS:
      i = sscanf(buf," %d %d %20s %f %f",
                 &nbfunc,&combrule,genpairs,&fudgeLJ,&fudgeQQ);
      if(i < 3) { // didn't get enough parameters 
        fprintf(stderr,"syntax error in DEFAULTS\n");
        exit(1);
      }
      if(nbfunc != 1) { // I don't know how it works with nbfunc=2
        fprintf(stderr,"Non-bonded function != 1 unsupported in DEFAULTS\n");
        exit(1);
      }
      if(combrule != 1) { // same here
        fprintf(stderr,"Combination rule != 1 unsupported in DEFAULTS\n");
        exit(1);
      }
      if(0==strcmp(genpairs,"yes")) {
        genPairs=1;
        if(i!=5) {
          fprintf(stderr,"syntax error in DEFAULTS\n");
          exit(1);
        }
        // else fudgeLJ and fudgeQQ got written automatically
      }
      else genPairs=0;

      break;

    case NONBOND:
      if(5 != sscanf(buf," %5s %5s %d %f %f",
                     typea, typeb, &funct, &c6, &c12)) {
        fprintf(stderr,"Syntax error in NONBOND\n");
        exit(1);
      }
      // convert kJ/mol*nm6 to kcal/mol*A6 and ..12 ..12
      c6 =  c6/JOULES_PER_CALORIE*1E6;
      c12= c12/JOULES_PER_CALORIE*1E12;
      vdwTable.addType(typea,typeb,c6,c12);
      break;

    case BONDS:
      i = sscanf(buf," %d %d %d %f %f",
                 &atomi,&atomj,&funct,&c0,&c1);
      atomi--; // shift right away to a zero-indexing
      atomj--;
      if(i==3) {
       tmptypea = genericMols[genericMols.size()-1]->getAtom(atomi)->getType();
       tmptypeb = genericMols[genericMols.size()-1]->getAtom(atomj)->getType();
        /* find the index and parameters */
        index = bondTable.getParams(tmptypea, tmptypeb, funct, &b0, &kB);
        if(index==-1) {
          fprintf(stderr,"Required bondtype %s--%s (function %d) not found.\n",
                  tmptypea,tmptypeb,funct);
          exit(1);
        }
      }
      else if(i==5) {
        /* first set the values of b0 and kB correctly */
        b0 = c0*10; /* convert nm to A */
        if(funct==1) { /* harmonic potential */
          /* convert kJ/nm2 to kcal/A2 and use E=kx2 instead of half that. */
          kB = c1/JOULES_PER_CALORIE/100/2;
        }
        else if(funct==2) { /* special fourth-order potential */
          /* convert to the normal harmonic constant and kJ/nm2 to kcal/A2 */
          kB = 2*c1*c0*c0/JOULES_PER_CALORIE/100;
          kB /= 2; /* use the NAMD system where E=kx2 */
        }
        else {
          fprintf(stderr,"I don't know what funct=%d means in BONDS\n",funct);
          exit(1);
        }
        /* look up the index */
        index = bondTable.getIndex(b0,kB,funct);
      }
      else {
        fprintf(stderr,"Syntax error in BONDS\n");
        exit(1);
      }

      genericMols[genericMols.size()-1]->addBond(atomi,atomj,index);
      break;
      
    case BONDTYPES:
      if(5 != sscanf(buf," %5s %5s %d %f %f",
                     typea,typeb,&funct,&c0,&c1)) {
        fprintf(stderr,"Syntax error in BONDTYPES\n");
        exit(1);
      }

      /* first set the values of b0 and kB correctly */
      b0 = c0*10; /* convert nm to A */
      if(funct==1) { /* harmonic potential */
        /* convert kJ/nm2 to kcal/A2 and use E=kx2 instead of half that. */
        kB = c1/JOULES_PER_CALORIE/100/2;
      }
      else if(funct==2) { /* special fourth-order potential */
        /* convert to the normal harmonic constant and kJ/nm2 to kcal/A2 */
        kB = 2*c1*c0*c0/JOULES_PER_CALORIE/100;
        kB /= 2; /* use the NAMD system where E=kx2 */
      }
      else {
 fprintf(stderr,"I don't know what funct=%d means in BONDTYPES\n",funct);
        exit(1);
      }

      bondTable.addType(typea,typeb,b0,kB,funct);
      break;

    case ANGLES:
      i = sscanf(buf," %d %d %d %d %f %f",
                     &atomi,&atomj,&atomk,&funct,&c0,&c1);
      atomi--; // shift right away to a zero-indexing
      atomj--;
      atomk--;
      if(i == 4) { /* we have to look up the last two parameters */
       tmptypea = genericMols[genericMols.size()-1]->getAtom(atomi)->getType();
       tmptypeb = genericMols[genericMols.size()-1]->getAtom(atomj)->getType();
       tmptypec = genericMols[genericMols.size()-1]->getAtom(atomk)->getType();
        /* find the index and parameters */
        index = angleTable.getParams(tmptypea, tmptypeb, tmptypec,
                                        funct, &th0, &kth);
        if(index==-1) {
          fprintf(stderr,
                  "Required angletype %s--%s--%s (function %d) not found.\n",
                  tmptypea,tmptypeb,tmptypec,funct);
          exit(1);
        }
      }
      else if(i == 6) {
        /* first set the values of th0 and kth correctly */
        if(funct == 1) {
          th0 = c0; /* both are in degrees */
          kth = c1/JOULES_PER_CALORIE/2; /* convert kJ/rad2 to kcal/rad2 and use E=kx2 */
        }
        else if(funct == 2) {
          th0 = c0; /* both are in degrees */
          /* convert G96 kJ to kcal/rad2 and use E=kx2 */
          kth = sin(th0*PI/180)*sin(th0*PI/180)*c1/JOULES_PER_CALORIE/2;
        }
        else {
          fprintf(stderr,"I don't know what funct=%d means in ANGLES\n",funct);
          exit(1);
        }
        /* add the angle type to our table */
        index = angleTable.getIndex(th0,kth,funct);
      }
      else {
        fprintf(stderr,"Syntax error (%d args) in ANGLES: %s\n",i,buf);
        exit(1);
      }

      /* add the angle to our table */
      genericMols[genericMols.size()-1]->addAngle(atomi,atomj,atomk,index);
      break;

    case ANGLETYPES:
      if(6 != sscanf(buf," %5s %5s %5s %d %f %f",
                     typea,typeb,typec,&funct,&c0,&c1)) {
        fprintf(stderr,"Syntax error in ANGLETYPES\n");
        exit(1);
      }
      /* first set the values of th0 and kth correctly */
      if(funct == 1) {
        th0 = c0; /* both are in degrees */
        kth = c1/JOULES_PER_CALORIE/2; /* convert kJ/rad2 to kcal/rad2 and use E=kx2 */
      }
      else if(funct == 2) {
        th0 = c0; /* both are in degrees */
        /* convert G96 kJ to kcal/rad2 and use E=kx2 */
        kth = sin(th0*PI/180)*sin(th0*PI/180)*c1/JOULES_PER_CALORIE/2;
      }
      else {
        fprintf(stderr,"I don't know what funct=%d means in ANGLETYPES\n",
                funct);
        exit(1);
      }
      angleTable.addType(typea,typeb,typec,th0,kth,funct);
      break;

    case DIHEDRALS:
      i = sscanf(buf," %d %d %d %d %d %f %f %f %f %f %f",
                 &atomi,&atomj,&atomk,&atoml,&funct,
                 &c[0],&c[1],&c[2],&c[3],&c[4],&c[5]);
      atomi--; // shift right away to a zero-indexing
      atomj--;
      atomk--;
      atoml--;
      if(i==5) { /* we have to look up the parameters */
       tmptypea = genericMols[genericMols.size()-1]->getAtom(atomi)->getType();
       tmptypeb = genericMols[genericMols.size()-1]->getAtom(atomj)->getType();
       tmptypec = genericMols[genericMols.size()-1]->getAtom(atomk)->getType();
       tmptyped = genericMols[genericMols.size()-1]->getAtom(atoml)->getType();
        /* find the index and parameters */
        index = dihedralTable.getParams(tmptypea, tmptypeb, tmptypec,
                                        tmptyped, funct, c, &mult);
        if(index==-1) {
          fprintf(stderr,
             "Required dihedraltype %s--%s--%s--%s (function %d) not found.\n",
                  tmptypea,tmptypeb,tmptypec,tmptyped,funct);
          exit(1);
        }
      }
      else if(i==7 || i==8 || i==11) { /* the parameters are given */
        if(funct==1 || funct==2) { /* we should have two parameters */
          if(i!=7+(funct==1)) {    /* plus a multiplicity for funct==1 */
            fprintf(stderr,"Must have 7 args for funct=1,2\n");
            exit(1);
          }
          c[0] = c[0]; /* both in deg */
          c[1] = c[1]/JOULES_PER_CALORIE; /* convert kJ to kcal and still use E=kx2*/
          /* for funct==1 these are both divided by rad^2 */
          if(funct==1) {
            mult=(int)(c[2]+0.5); /* round to nearest integer :p */
          }
        }
        else if(funct==3) {
          if(i!=11){
            fprintf(stderr,"Must have 11 args for funct=3\n");
            exit(1);
          }

          for(j=0;j<6;j++) {
            c[j] = c[j]/JOULES_PER_CALORIE; /* convert kJ to kcal and E=Cn cos^n */
          }
        }
        else {
          fprintf(stderr,
                  "I don't know what funct=%d means in DIHEDRALS\n",funct);
          exit(1);
        }
        index = dihedralTable.getIndex(c,mult,funct);
      }
      else {
        fprintf(stderr,"Syntax error (%d args) in DIHEDRALS: %s\n",i,buf);
        exit(1);
      }

      /* add the dihedrals to our table */
      genericMols[genericMols.size()-1]->addDihedral(atomi,atomj,atomk,atoml,
                                                     index);
      break;
    case DIHEDRALTYPES:
      i = sscanf(buf," %5s %5s %d %f %f %f %f %f %f",
                 typea,typeb,&funct,&c[0],&c[1],&c[2],&c[3],&c[4],&c[5]);
      if(funct == 1 || funct == 2) {
        if(i!=5+(funct==1)) { /* 6 for f=2, 5 for f=1 */
          fprintf(stderr,"Syntax error in DIHEDRALTYPES: %s\n",buf);
          exit(1);
        }
        c[0] = c[0]; /* both in deg */
        c[1] = c[1]/JOULES_PER_CALORIE; /* convert kJ to kcal and still use E=kx2*/
        /* for funct==1 these are both divided by rad^2 */
        if(funct==1) {
          mult=(int)(c[2]+0.5); /* round to nearest integer :p */
        }
      }
      else if(funct == 3) {
        if(i!=9) {
          fprintf(stderr,"Syntax error in DIHEDRALTYPES\n");
          exit(1);
        }
        for(j=0;j<6;j++) {
          c[j] = c[j]/JOULES_PER_CALORIE; /* convert kJ to kcal and E=Cn cos^n */
        }
      }
      else {
        fprintf(stderr,"I don't know what funct=%d means in DIHEDRALTYPES\n",
                funct);
        exit(1);
      }
      dihedralTable.addType(typea,typeb,c,mult,funct);
      break;
      
    case ATOMS:
      i = sscanf(buf," %d %5s %d %5s %5s %d %f %f",
                   &num, type, &resnum, restype,
                 atomname, &chargegp, &charge, &mass);
      if(i==7) { /* XXX temporary - I should be able to get more
                    params */
        typenum = atomTable.getParams(type,&mass,&junkf,&junkf,&junkf);
        i=8;
      }
      else {
        if(i!=8) {
          fprintf(stderr,"Syntax error in ATOMS\n");
          exit(1);
        }
        // just get the type number
        typenum = atomTable.getParams(type,&junkf,&junkf,&junkf,&junkf);
      }
      genericMols[genericMols.size()-1]->addAtom(type,typenum,resnum,
                                                 restype,atomname,charge,mass);
      break;

    case ATOMTYPES:
      if(6 != sscanf(buf," %5s %f %f %5s %f %f",type,&mass,&charge,
                     particletype,&c6,&c12)) {
        fprintf(stderr,"Syntax error in ATOMTYPES: %s\n",buf);
        exit(1);
      }
      /* conversions:
         c6  - kJ/mol nm6  -> kcal/mol A6
         c12 - kJ/mol nm12 -> kcal/mol A12 */
      atomTable.addType(type,mass,charge,
                            c6/JOULES_PER_CALORIE*1E6,
                            c12/JOULES_PER_CALORIE*1E12);
      break;

    case MOLECULETYPE:
      if(2!=sscanf(buf," %20s %d",molname,&nrexcl)) {
        fprintf(stderr,"Syntax error in MOLECULETYPE: %s\n",buf);
        exit(1);
      }

      /* add another generic molecule holder */
      genericMols.add(new GenericMol(molname));
      break;

    case MOLECULES:
      if(2!=sscanf(buf," %20s %d",molname,&copies)) {
        fprintf(stderr,"Syntax error in MOLECULES: %s\n",buf);
        exit(1);
      }
      
      /* search for the specified molecule and make a molInst of it*/
      moleculeinstance = NULL;
      for(i=0;i<genericMols.size();i++) {
        if(0==strcmp(molname,genericMols[i]->getName())) {
          moleculeinstance = new MolInst(genericMols[i],copies);
          break;
        }
      }

      if(moleculeinstance==NULL) {
        fprintf(stderr,"Molecule %s undefined.\n",molname);
        exit(1);
      }
      
      /* put it at the end of the list */
      molInsts.add(moleculeinstance);

      break;
    }
  }

  fclose(f);
}

/* returns the number of bonds in the file */
int GromacsTopFile::getNumBonds() const {
  int n=0,i;
  for(i=0;i<molInsts.size();i++) {
    n += molInsts[i]->getNum() *
         molInsts[i]->getMol()->getNumBonds();
  }
  return n;
}

/* returns the number of angles in the file */
int GromacsTopFile::getNumAngles() const {
  int n=0,i;
  for(i=0;i<molInsts.size();i++) {
    n += molInsts[i]->getNum() *
         molInsts[i]->getMol()->getNumAngles();
  }
  return n;
}

/* returns the number of dihedral angles in the file */
int GromacsTopFile::getNumDihedrals() const {
  int n=0,i;
  for(i=0;i<molInsts.size();i++) {
    n += molInsts[i]->getNum() *
         molInsts[i]->getMol()->getNumDihedrals();
  }
  return n;
}

/* getBond puts the information about bond number <num> into the
   spaces pointed to by the other arguments.  Bond number 0 is the
   first bond in the list.
   For atomi and atomj, 1 is the first atom in the list. */
void GromacsTopFile::getBond(int n, int *atomi, int *atomj, int *bondtype) const {
  /* figure out what instance this bond is in */
  int nbonds=0;
  int natoms=0;
  int i;
  for(i=0;i<molInsts.size();i++) {
    int molbonds = molInsts[i]->getMol()->getNumBonds();
    int molatoms = molInsts[i]->getMol()->getNumAtoms();
    int newbonds = molInsts[i]->getNumBonds();
    int newatoms = molInsts[i]->getNumAtoms();
    
    if(nbonds+newbonds-1 >= n) {
      /* the bond is in this MolInst */
      int localnum = (n-nbonds) % molbonds;   /* number within this inst */
      int instnum = (n-nbonds) / molbonds;    /* number of instances before */
      int addatoms = natoms+instnum*molatoms; /* extra atoms to add to atomi */

      const GenericBond *b = molInsts[i]->getMol()->getBond(localnum);

      *atomi = b->getAtomi() + addatoms;
      *atomj = b->getAtomj() + addatoms;
      *bondtype = b->getType();
      break;
    }

    nbonds += newbonds;
    natoms += newatoms;
  }
}

/* getAngle puts the information about angle number <num> into the
   spaces pointed to by the other arguments.  Angle number 0 is the
   first angle in the list.
*/
void GromacsTopFile::getAngle(int n, int *atomi, int *atomj, int *atomk,
                          int *angletype) const {
  /* figure out what instance this angle is in */
  int nangles=0;
  int natoms=0;
  int i;
  for(i=0;i<molInsts.size();i++) {
    int molangles = molInsts[i]->getMol()->getNumAngles();
    int molatoms = molInsts[i]->getMol()->getNumAtoms();
    int newangles = molInsts[i]->getNumAngles();
    int newatoms = molInsts[i]->getNumAtoms();
    
    if(nangles+newangles-1 >= n) {
      /* the angle is in this MolInst */
      int localnum = (n-nangles) % molangles; /* number within this inst */
      int instnum = (n-nangles) / molangles;  /* number of instances before */
      int addatoms = natoms+instnum*molatoms; /* extra atms to add to atmi */

      const GenericAngle *a = molInsts[i]->getMol()->getAngle(localnum);

      *atomi = a->getAtomi() + addatoms;
      *atomj = a->getAtomj() + addatoms;
      *atomk = a->getAtomk() + addatoms;
      *angletype = a->getType();
      break;
    }

    nangles += newangles;
    natoms += newatoms;
  }
}

/* getDihedral puts the information about dihedral number <num> into
   the spaces pointed to by the other arguments.  Dihedral number 0
   is the first angle in the list. */
void GromacsTopFile::getDihedral(int n, int *atomi, int *atomj, int *atomk,
                              int *atoml, int *type) const {
  /* figure out what instance this angle is in */
  int ndihedrals=0;
  int natoms=0;
  int i;
  for(i=0;i<molInsts.size();i++) {
    int moldihedrals = molInsts[i]->getMol()->getNumDihedrals();
    int molatoms = molInsts[i]->getMol()->getNumAtoms();
    int newdihedrals = molInsts[i]->getNumDihedrals();
    int newatoms = molInsts[i]->getNumAtoms();
    
    if(ndihedrals+newdihedrals-1 >= n) {
      /* the dihedral is in this MolInst */
      int localnum = (n-ndihedrals) % moldihedrals; /* number in this inst */
      int instnum = (n-ndihedrals) / moldihedrals; /* number of insts before */
      int addatoms = natoms+instnum*molatoms; /*extra atms to add to atmi*/

      const GenericDihedral *a = molInsts[i]->getMol()->getDihedral(localnum);

      *atomi = a->getAtomi() + addatoms;
      *atomj = a->getAtomj() + addatoms;
      *atomk = a->getAtomk() + addatoms;
      *atoml = a->getAtoml() + addatoms;
      *type = a->getType();
      break;
    }

    ndihedrals += newdihedrals;
    natoms += newatoms;
  }
}

/* getNumAtoms returns the number of atoms stored in the file. */
int GromacsTopFile::getNumAtoms() const {
  int n=0;
  int i;
  for(i=0;i<molInsts.size();i++) {
    n += molInsts[i]->getNum() *
         molInsts[i]->getMol()->getNumAtoms();
  }
  return n;
}

/* getAtom puts the information about atom number <n> into the
   spaces pointed to by the other arguments.  The string buffers
   must be at least 11 characters long. */
void GromacsTopFile::getAtom(int num,
               int *residue_number, char *residue_name,
               char *atom_name, char *atom_type, int *atom_typenum,
               Real *charge, Real *mass) 
  const {
  int natoms=0,n=num; // zero-indexed array
  int resnum=0;
  /* figure out what instance the atom is in, and what residue number
     it has */
  int i;

  for(i=0;i<molInsts.size();i++) {
    int numnew = molInsts[i]->getNumAtoms(); /* atoms in this MolInst */
    int resmol = molInsts[i]->getMol()->getNumRes(); /* # res/mol */
    int newres = molInsts[i]->getNumRes(); /* residues in this MolInst */

    if(natoms+numnew-1 >= n) {
      /* the atom is in this molInst */
      int localnum = (n-natoms) % molInsts[i]->getMol()->getNumAtoms();
      int instnum  = (n-natoms) / molInsts[i]->getMol()->getNumAtoms();
      
      // getAtom is zero-indexed
      const GenericAtom *a = molInsts[i]->getMol()->getAtom(localnum);
      int residue = resnum + resmol*instnum + a->getResNum();

      *residue_number = residue;
      strncpy(residue_name,a->getResType(),11);
      strncpy(atom_name,a->getAtomName(),11);
      strncpy(atom_type,a->getType(),11);
      *charge=a->getCharge();
      *mass=a->getMass();
      *atom_typenum=a->getTypeNum();
      break;
    }

    /* go to the next molInst */
    natoms += numnew;
    resnum += newres;
  }
}

GenericAtom::GenericAtom(const char *theType, int theTypeNum, int theResNum,
                         const char *theResType,const char *theAtomName,
                         Real theCharge, Real theMass) {
  strncpy(type,theType,NAMESIZE+1);
  typeNum = theTypeNum;
  resNum = theResNum;
  strncpy(resType,theResType,NAMESIZE+1);
  strncpy(atomName,theAtomName,NAMESIZE+1);
  charge = theCharge;
  mass = theMass;
}

/* initializes this to be a bond between atoms <i> and <j> of type
   <type>, with <next> pointing to the next bond in the list. */
GenericBond::GenericBond(int i, int j, int theType) {
  atomi=i;
  atomj=j;
  type=theType;
}

/* initializes this to be a angle between atoms <i>, <j>, and <k> of
   type <type>, with <next> pointing to the next angle in the list. */
GenericAngle::GenericAngle(int i, int j, int k, int theType) {
  atomi=i;
  atomj=j;
  atomk=k;
  type=theType;
}

/* initializes this to be a angle between atoms <i>, <j>, <k>, and
   <l> of type <type> */
GenericDihedral::GenericDihedral(int i, int j, int k, int l, int theType) {
  atomi=i;
  atomj=j;
  atomk=k;
  atoml=l;
  type=theType;
}

/* adds a bond to the list */
void GenericMol::addBond(int atomi, int atomj, int type) {
  bondList.add(new GenericBond(atomi, atomj, type));
}

/* adds an angle to the list */
void GenericMol::addAngle(int atomi, int atomj, int atomk, int type) {
  angleList.add(new GenericAngle(atomi, atomj, atomk, type));
}

/* adds a dihedral to the list */
void GenericMol::addDihedral(int atomi, int atomj, int atomk,
                             int atoml, int type) {
  dihedralList.add(new GenericDihedral(atomi, atomj, atomk,
                                             atoml, type));
}

/* adds an atom to the list */
void GenericMol::addAtom(const char *theType, int theTypeNum, int theResNum,
                         const char *theResType,const char *theAtomName,
                         Real theCharge, Real theMass) {

  atomList.add(new GenericAtom(theType,theTypeNum,theResNum,theResType,
                           theAtomName,theCharge,theMass));
}

/* initializes this to be the molecule with name <name> */
GenericMol::GenericMol(const char *theName) {
  name = strdup(theName);
}

/* gets a bond */
const GenericBond *GenericMol::getBond(int n) const {
  /* double-check */
  if(n >= bondList.size() || n<0) {
    fprintf(stderr,"Bond index %d out of bounds.\n",n);
    exit(1);
  }
  return bondList[n];
}

/* gets a angle */
const GenericAngle *GenericMol::getAngle(int n) const {
  /* double-check */
  if(n >= angleList.size() || n<0) {
    fprintf(stderr,"Angle index %d out of bounds.\n",n);
    exit(1);
  }
  return angleList[n];
}

/* gets a dihedral */
const GenericDihedral *GenericMol::getDihedral(int n) const {
  /* double-check */
  if(n >= dihedralList.size() || n<0) {
    fprintf(stderr,"Dihedral index %d out of bounds.\n",n);
    exit(1);
  }
  return dihedralList[n];
}

/* gets an atom */
const GenericAtom *GenericMol::getAtom(int n) const {
  /* double-check */
  if(n >= atomList.size() || n<0) {
    fprintf(stderr,"Atom index %d out of bounds for %s.\n",n,name);
    exit(1);
  }
  return atomList[n];
}

/* initializes this to represent <theNum> copies of <theMol> */
MolInst::MolInst(const GenericMol *theMol, int theNum) {
  mol = theMol;
  num = theNum;
}

/* get the total number of various things here */
int MolInst::getNumAtoms() const {
  return getMol()->getNumAtoms() * getNum();
}
int MolInst::getNumBonds() const {
  return getMol()->getNumBonds() * getNum();
}
int MolInst::getNumAngles() const {
  return getMol()->getNumAngles() * getNum();
}
int MolInst::getNumDihedrals() const {
  return getMol()->getNumDihedrals() * getNum();
}
int MolInst::getNumRes() const {
  return getMol()->getNumRes()
    * getNum();
}

/* this gets the index of a bond in the table (adding an entry if
   none exists).
   b0: the natural bond length in A.
   kB: the spring constant in kcal/mol/A^2, where E=kx^2 to 1st order.
   funct: 1 for harmonic potential, 2 for fourth-order GROMACS96
   potential. */
int BondTable::getIndex(float b0, float kB, int funct) {
  /* check to see if it is in the table already */
  int i;
  for(i=0;i<b0Array.size();i++) {
    if(fabs(b0-b0Array[i])<0.00001 &&
       fabs(kB-kBArray[i])<0.00001 &&
       funct == functArray[i]) {
      return i;
    }
  }
  
  /* nope, it wasn't in the table add a new element! */
  b0Array.add(b0);
  kBArray.add(kB);
  functArray.add(funct);
  typeaArray.add(NULL);
  typebArray.add(NULL);
  return b0Array.size()-1;
}

/* this gets the index of a angle in the table (adding an entry if
   none exists).
   b0: the natural angle length in A.
   kB: the spring constant in kcal/mol/A^2, where E=kx^2 to 1st order.
   funct: 1 for harmonic potential, 2 for fourth-order GROMACS96
   potential. */
int AngleTable::getIndex(float th0, float kth, int funct) {
  /* check to see if it is in the table already */
  int i;
  for(i=0;i<th0Array.size();i++) {
    if(fabs(th0-th0Array[i])<0.00001 &&
       fabs(kth-kthArray[i])<0.00001 &&
       funct == functArray[i]) {
      return i;
    }
  }
  
  /* nope, it wasn't in the table add a new element! */
  th0Array.add(th0);
  kthArray.add(kth);
  functArray.add(funct);
  typeaArray.add(NULL);
  typebArray.add(NULL);
  typecArray.add(NULL);
  return th0Array.size()-1;
}

/* this function gets the index of a dihedral angle in the table
   (adding an entry if none exists).  The required number of
   parameters (see notes on class DihedralTable) must be stored in
   the array <c> */
int DihedralTable::getIndex(const float *c, int mult, int funct) {
  /* check to see if it is in the table already */
  int i,j,jmax;
  if(funct==1 || funct==2) { /* for these we only need two params */
    jmax=2;
  } else { /* for RB we need all 6 */
    jmax=6;
  }

  for(i=0;i<cArray.size();i++) {
    int mismatch=0;
    if(mult != multArray[i]) continue;
    if(funct != functArray[i]) continue;

    for(j=0;j<jmax;j++) {
      if(fabs(c[j]-cArray[i][j])>=0.00001) {
        mismatch=1;
        break;
      }
    }
    if(!mismatch) {
      /* all of the parameters matched */
      return i;
    }
  }
  
  /* nope, it wasn't in the table add a new element! */
  addType(NULL,NULL,c,mult,funct);
  return cArray.size()-1;
}

/* getBondParams puts the parameters for bond-type <num> into the
   spaces pointed to by the other arguments. 
   b0 - natural length in A
   kB - spring constant in kcal/A^2 - E=kx^2
   funct - 1 for harmonic, 2 for special fourth-order GROMOS96 thing */
void GromacsTopFile::getBondParams(int num, float *b0, float *kB, int *funct)
  const {
  bondTable.getParams(num,b0,kB,funct);
}

/* getAngleParams puts the parameters for angle-type <num> into the
   spaces pointed to by the other arguments.  */
void GromacsTopFile::getAngleParams(int num, float *th0, float *kth, int *funct)
  const {
  angleTable.getParams(num,th0,kth,funct);
}

void GromacsTopFile::getDihedralParams(int num, float *c, int *mult, int *funct)
  const {
  dihedralTable.getParams(num,c,mult,funct);
}

/* getParams puts the parameters for bond-type <num> into the
   spaces pointed to by the other arguments. 
   b0 - natural length in A
   kB - spring constant in kcal/A^2 - E=kx^2
   funct - 1 for harmonic, 2 for special fourth-order GROMOS96 thing */
void BondTable::getParams(int num, float *b0, float *kB, int *funct)
  const {
  *b0=b0Array[num];
  *kB=kBArray[num];
  *funct=functArray[num];
}

/* getParams puts the parameters for angle-type <num> into the
   spaces pointed to by the other arguments.  */
void AngleTable::getParams(int num, float *th0, float *kth, int *funct)
  const {
  *th0=th0Array[num];
  *kth=kthArray[num];
  *funct=functArray[num];
}

/* getParams puts the parameters for angle-type <num> into the
   spaces pointed to by the other arguments.  The required number of
   parameters (see notes on class DihedralTable) will be stored in
   the array <c>, so make sure that c has size >= 6 */
void DihedralTable::getParams(int num, float *c, int *mult, int *funct) const {
  int i;

  *funct=functArray[num]; /* first set the function */

  if(*funct==1 || *funct==2) { /* for these we only need two params */
    c[0]=cArray[num][0];
    c[1]=cArray[num][1];
  } else if(*funct==3) { /* for RB we need all 6 */
    for(i=0;i<6;i++) {
      c[i]=cArray[num][i];
    }
  } else { /* error */
    fprintf(stderr,"Bad function number %d - don't know what to do!\n",
            funct);
    exit(1);
  }

  if(*funct==1) { /* return the multiplicity */
    *mult=multArray[num];
  }
}

/* this adds a entry for a angle type to the table, including the
   three atoms involved in the angle */
void AngleTable::addType(const char *typea, const char *typeb,
                  const char *typec, float th0, float kth, int funct) {
  typeaArray.add(strdup(typea));
  typebArray.add(strdup(typeb));
  typecArray.add(strdup(typec));
  th0Array.add(th0);
  kthArray.add(kth);
  functArray.add(funct);
}

/* this adds a entry for a angle type to the table, including the
   two atoms involved in the angle.  The required number of
   parameters (see notes on class DihedralTable) must be stored in
   the array <c>.  Note that these two angles really refer to either
   atoms A and D or B and C depending on the dihedral type.  */
void DihedralTable::addType(const char *typea, const char *typeb,
                         const float *c, int mult, int funct) {
  float *cNew;
  int i;

  if(typea != NULL) typeaArray.add(strdup(typea));
  if(typeb != NULL) typebArray.add(strdup(typeb));
  functArray.add(funct);

  if(funct==1) multArray.add(mult);
  else multArray.add(0);

  if(funct==1 || funct==2) { /* for these we only need two params */
    cNew = new float[2];
    cNew[0]=c[0];
    cNew[1]=c[1];
  } else if(funct==3) { /* for RB we need all 6 */
    cNew = new float[6];
    for(i=0;i<6;i++) {
      cNew[i]=c[i];
    }
  } else { /* error */
    fprintf(stderr,"Bad function number %d - don't know what to do!\n",
            funct);
    exit(1);
  }
  
  /* now push the actual parameters */
  cArray.add(cNew);
}

/* This version looks up the angle by atom type - the direction of
   the types doesn't matter (it can be A--B--C or C--B--A)!  If the
   specified angle/function combination cannot be found, it returns
   -1, otherwise it returns the index of the angletype. */
int AngleTable::getParams(const char *typea, const char *typeb,
                   const char *typec, int funct,
                   float *th0, float *kth) const {
  int i;
  for(i=0;i<th0Array.size();i++) {
    if(typeaArray[i] == NULL || typebArray[i] == NULL
       || typecArray[i] == NULL) continue;
    if( (0==strcmp(typea,typeaArray[i]) &&  /* A--B--C */
         0==strcmp(typeb,typebArray[i]) &&
         0==strcmp(typec,typecArray[i]))    /* or */
     || (0==strcmp(typec,typeaArray[i]) &&
         0==strcmp(typeb,typebArray[i]) &&  /* C--B--A */
         0==strcmp(typea,typecArray[i]))
     && funct==functArray[i]) {
      *th0 = th0Array[i];
      *kth = kthArray[i];
      return i;
    }
  }
  return -1;
}

/* see the (long) notes in the prototype definition */
int DihedralTable::getParams(const char *typea, const char *typeb,
                             const char *typec, const char *typed,
                             int funct, float *c, int *mult) const {
  int i,j;
  const char *comparea, *compareb;
  
  if(funct == 1 || funct == 3) { /* for these, we use the inner atoms */
    comparea = typeb;
    compareb = typec;
  } else { /* use the outer atoms */
    comparea = typea;
    compareb = typed;
  }

  for(i=0;i<cArray.size();i++) {
    if(typeaArray[i] == NULL || typebArray[i] == NULL)
      continue; /* no atom types specified */
    if(functArray[i] != funct)
      continue; /* wrong function type */

    if( (0==strcmp(comparea,typeaArray[i]) &&  /* A--B */
         0==strcmp(compareb,typebArray[i]))    /*  or  */
     || (0==strcmp(compareb,typeaArray[i]) &&
         0==strcmp(comparea,typebArray[i]))    /* B--A */
          ) {
      if(funct==1 || funct==2) { /* for these we only need two params */
        c[0]=cArray[i][0];
        c[1]=cArray[i][1];
        if(funct==1) {
          *mult = multArray[i];
        }
      } else if(funct==3) { /* for RB we need all 6 */
        for(j=0;j<6;j++) {
          c[j]=cArray[i][j];
        }
      }
      return i;
    }
  }
  return -1;
}

/* this adds a entry for a bond type to the table, including the two
   atoms involved in the bond */
void BondTable::addType(const char *typea, const char *typeb,
                            float b0, float kB, int funct) {
  b0Array.add(b0);
  kBArray.add(kB);
  functArray.add(funct);
  typeaArray.add(strdup(typea));
  typebArray.add(strdup(typeb));
}

/* This version looks up the bond by atom type - the order of the
   types doesn't matter! */
int BondTable::getParams(const char *typea, const char *typeb,
                              int funct, float *b0, float *kB) const {
  int i;
  for(i=0;i<b0Array.size();i++) {
    if(typeaArray[i] == NULL || typebArray[i] == NULL) continue;
    if( (0==strcmp(typea,typeaArray[i]) &&
         0==strcmp(typeb,typebArray[i]))
     || (0==strcmp(typeb,typeaArray[i]) &&
         0==strcmp(typea,typebArray[i]))
     && funct==functArray[i]) {
      *b0 = b0Array[i];
      *kB = kBArray[i];
      return i;
    }
  }
  return -1;
}

/* this adds a entry for an atom type to the table */
void AtomTable::addType(const char *type, float m, float q,
                            float c6, float c12) {
  typeArray.add(strdup(type));
  mArray.add(m);
  qArray.add(q);
  c6Array.add(c6);
  c12Array.add(c12);
}

/* This looks up the atom type by number, returning it in the string
   <type>, which must be at least 11 characters long. */
void AtomTable::getType(int num, char *type) const {
  if(num>=mArray.size() || num<0) {
    fprintf(stderr,"atomParam index %d out of bounds!\n",num);
    exit(1);
  }
  strncpy(type,typeArray[num],NAMESIZE+1);
}

/* This looks up the atom by type - if it cannot be found, this
   returns -1, otherwise this returns the index of the atomtype (for
   consistency - this is really a useless number.) */
int AtomTable::getParams(const char *type, float *m, float *q,
                  float *c6, float *c12) const {
  int i;
  for(i=0;i<mArray.size();i++) {
    if(typeArray[i]==NULL) {
      fprintf(stderr,"Found NULL atom type in list.\n");
      exit(1);
    }
    if(0==strcmp(typeArray[i],type)) {
      *m = mArray[i];
      *q = qArray[i];
      *c6 = c6Array[i];
      *c12 = c12Array[i];
      return i;
    }
  }
  return -1;
}

/* finds the index from the two interacting types, or returns -1 */
int VDWTable::getIndex(const char *typea, const char *typeb) const {
  int i;
  for(i=0;i<c6Array.size();i++) {
    if((0==strcmp(typea,typeAArray[i]) &&
        0==strcmp(typeb,typeBArray[i])) ||
       (0==strcmp(typeb,typeAArray[i]) &&
        0==strcmp(typea,typeBArray[i]))) {
      return i;
    }
  }
  return -1;
}

/* these add VDW parameters to the list */
void VDWTable::addType(const char *typea, const char *typeb, Real c6, 
                       Real c12) {
  int i;

  /* check to see if the pair is already in the table */
  i = getIndex(typea,typeb);
  if(i != -1) {  /* it was in the table */
    c6Array[i] = c6;
    c12Array[i] = c12;
  }
  else { /* it wasn't in the table - add it! */
    typeAArray.add(strdup(typea));
    typeBArray.add(strdup(typeb));
    c6Array.add(c6);
    c12Array.add(c12);
    c6PairArray.add(0);
    c12PairArray.add(0);
  }
}

void VDWTable::add14Type(const char *typea, const char *typeb,
               Real c6pair, Real c12pair) {
  int i;

  /* check to see if the pair is already in the table */
  i = getIndex(typea,typeb);
  if(i != -1) {  /* it was in the table */
    c6PairArray[i] = c6pair;
    c12PairArray[i] = c12pair;
  }
  else { /* it wasn't in the table - add it! */
    typeAArray.add(strdup(typea));
    typeBArray.add(strdup(typeb));
    c6PairArray.add(c6pair);
    c12PairArray.add(c12pair);
    c6Array.add(0);
    c12Array.add(0);
  }
}

/* this returns the VDW interaction parameters that were added to
   the list earlier, and returns the index or the parameters (a
   useless number) or -1 if it can't find the specified types. */
int VDWTable::getParams(const char *typea, const char *typeb,
              Real *c6, Real *c12, Real *c6pair, Real *c12pair) const {
  int i;
  /* check to see if the pair is already in the table */
  i = getIndex(typea,typeb);
  if(i != -1) {  /* it was in the table - return the parameters  */
    *c6 = c6Array[i];
    *c12 = c12Array[i];
    *c6pair = c6PairArray[i];
    *c12pair = c12PairArray[i];
  }
  return i;
}

/* getVDWParams returs the Lennard-Jones bonding parameters for the
   specified two atom types, and the modified bonding parameters
   for 1-4 L-J interactons (c6pair, c12pair). */
void GromacsTopFile::getVDWParams(int numa, int numb,
                  Real *c6, Real *c12, Real *c6pair, Real *c12pair) const {
  int i,ret;
  Real c6a,c6b,c12a,c12b, mjunk, qjunk;
  char typea[11]="",typeb[11]="";

  /* get the string names corresponding to numa and numb */
  getAtomParams(numa,typea);
  getAtomParams(numb,typeb);
  if(typea[0]==0) { /* found a bug in my use of strncpy here once */
    fprintf(stderr,"Failed to get name of atom %d\n",numa);
    exit(1);
  }
  if(typeb[0]==0) {
    fprintf(stderr,"Failed to get name of atom %d\n",numb);
    exit(1);
  }

  /* first try - use the VDW table */
  i = vdwTable.getParams(typea,typeb,c6,c12,c6pair,c12pair);

  if(i==-1) {
    /* fallback - use the individual atom's parameters */
    ret = atomTable.getParams(typea, &mjunk, &qjunk, &c6a, &c12a);
    if(ret==-1) {
      fprintf(stderr,"Couldn't find atom type %s\n",typea);
      exit(1);
    }
    ret = atomTable.getParams(typeb, &mjunk, &qjunk, &c6b, &c12b);
    if(ret==-1) {
      fprintf(stderr,"Couldn't find atom type %s\n",typeb);
      exit(1);
    }
    
    *c6  = (float)(sqrt(c6a * c6b));
    *c12 = (float)(sqrt(c12a*c12b));
  }

  /* XXX we need to put exclusions in here, including those implied by
     the use of RB. This may be a  little tricky because we will not
     be able to support any other kinds of dihedral bonds with the
     same end pair as any RB dihedral.  OR else, maybe NAMD will
     support the exclusions directly? */

  if(!genPairs) {
    return; /* we're done */
  }    

  *c6pair  = *c6  * fudgeLJ;
  *c12pair = *c12 * fudgeLJ;
  /* now we're really done */
}
    

---