#include <GromacsTopFile.h>

Public Member Functions
	GromacsTopFile (char *filename)

char *	getSystemName () const

int	getNumAtoms () const

int	getNumBonds () const

int	getNumAngles () const

int	getNumDihedrals () const

int	getNumAtomParams () const

int	getNumBondParams () const

int	getNumAngleParams () const

int	getNumDihedralParams () const

void	getAtom (int num, int residue_number, char residue_name, char atom_name, char atom_type, int atom_typenum, Real charge, Real *mass) const

void	getAtomParams (int num, char *type) const

int	getNumPair () const

int	getNumLJPair () const

int	getNumGaussPair () const

int	getNumExclusions () const

void	getExclusions (int , int ) const

void	getBond (int num, int atomi, int atomj, int *bondtype) const

void	getBondParams (int num, Real b0, Real kB, int *funct) const

void	getAngle (int num, int atomi, int atomj, int atomk, int angletype) const

void	getAngleParams (int num, Real th0, Real kth, int *funct) const

void	getDihedral (int num, int atomi, int atomj, int atomk, int atoml, int *type) const

void	getDihedralParams (int num, Real c, int mult, int *funct) const

void	getPairLJArrays2 (int indexA, int indexB, Real pairC6, Real pairC12)

void	getPairGaussArrays2 (int indexA, int indexB, Real gaussA, Real gaussMu1, Real gaussSigma1, Real gaussMu2, Real gaussSigma2, Real gaussRepulsive)

void	getVDWParams (int typea, int typeb, Real c6, Real c12, Real c6pair, Real c7) const

Detailed Description

Definition at line 479 of file GromacsTopFile.h.

Constructor & Destructor Documentation

◆ GromacsTopFile()

GromacsTopFile::GromacsTopFile ( char * filename )

Definition at line 65 of file GromacsTopFile.C.

References ResizeArray< Elem >::add(), PairTable::addPairGaussType2(), PairTable::addPairLJType2(), AtomTable::addType(), BondTable::addType(), AngleTable::addType(), DihedralTable::addType(), VDWTable::addType(), ANGLES, ANGLETYPES, ANGSTROMS_PER_NM, ATOMS, ATOMTYPES, BONDS, BONDTYPES, bool_negative_number_warning_flag, DEFAULTS, DIHEDRALS, DIHEDRALTYPES, EXCLUSIONS, exclusions_atom_i, exclusions_atom_j, BondTable::getIndex(), AngleTable::getIndex(), DihedralTable::getIndex(), AtomTable::getParams(), BondTable::getParams(), AngleTable::getParams(), DihedralTable::getParams(), iout, iWARN(), JOULES_PER_CALORIE, LINESIZE, LONGNAMESIZE, MOLECULES, MOLECULETYPE, NAMD_die(), NAMESIZE, NONBOND, numExclusion, PAIRS, PI, ResizeArray< Elem >::size(), SYSTEM, and UNKNOWN.

                                              {
   fudgeLJ = fudgeQQ = 1.0;
   /* 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;
       // JLai
       else if(0==strcmp(modename,"pairs")) mode = PAIRS;
       else if(0==strcmp(modename,"exclusions")) mode = EXCLUSIONS;
       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*ANGSTROMS_PER_NM; /* 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 */
           if(i==7) {
               c[1] = c[1]/(2*JOULES_PER_CALORIE); /* convert kJ to kcal and still use E=kx2*/
           } else if (i==8 || i==11) {
               c[1] = c[1]/(1*JOULES_PER_CALORIE); /* convert kJ to kcal and still use E=kx2*/
           } 
           //c[1] = c[1]/(2*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;
     case PAIRS:
       int indexA;
       int indexB;
       int pairFunction;
       Real fA;
       Real fB;
       Real fC;
       Real fD;
       Real fE;
       Real fF;
       
       int countVariables;
       countVariables = sscanf(buf," %d %d %d %f %f %f %f %f %f",&indexA,&indexB,&pairFunction,&fA,&fB,&fC,&fD,&fE,&fF);
 
       if ((countVariables >= 4 && countVariables >= 10)) {
         fprintf(stderr,"Syntax error in PAIRS: %s\n",buf);
         exit(1);
       }
       
       // Shift the atom indices to be zero-based
       indexA--;
       indexB--;
 
       // Make sure that the indexA is less than indexB
       /*if ( indexA > indexB ) {
           int tmpIndex = indexA;
           indexB = indexA;
           indexA = tmpIndex;
           }*/
 
       if (pairFunction == 1) {
         
         // LJ code
         fA = (fA/JOULES_PER_CALORIE)*1E6;
         fB= (fB/JOULES_PER_CALORIE)*1E12;
         pairTable.addPairLJType2(indexA,indexB,fA,fB);
       } else if (pairFunction == 5){
         
         // Bare Gaussian potential
         fA = (fA/JOULES_PER_CALORIE); //-->gaussA
         fB = (fB*ANGSTROMS_PER_NM); //-->gaussMu1
         if(fC == 0) {
           char buff[100];
           sprintf(buff,"GromacsTopFile.C::Attempting to load zero into gaussSigma.  Please check the pair: %s\n",buf);
           NAMD_die(buff);
         }
         if(fC < 0 && !bool_negative_number_warning_flag) {
           iout << iWARN << "Attempting to load a negative standard deviation into the gaussSigma.  Taking the absolute value of the standard deviation.";
           bool_negative_number_warning_flag = true;
         }
         fC = (fC*ANGSTROMS_PER_NM); //-->gaussSigma1
         fC = 1.0/(2 * fC * fC); // Normalizes sigma
         pairTable.addPairGaussType2(indexA,indexB,fA,fB,fC);
       } else if (pairFunction == 6) {
         
         // Combined Gaussian + repulsive r^-12 potential
         fA = (fA/JOULES_PER_CALORIE); //-->gaussA
         fB = (fB*ANGSTROMS_PER_NM); //-->gaussMu1
         if(fC == 0) {
           char buff[100];
           sprintf(buff,"GromacsTopFile.C::Attempting to load zero into gaussSigma.  Please check the pair: %s\n",buf);
           NAMD_die(buff);
         }
         if(fC < 0 && !bool_negative_number_warning_flag) {
           iout << iWARN << "Attempting to load a negative standard deviation into the gaussSigma.  Taking the absolute value of the standard deviation.";
           bool_negative_number_warning_flag = true;
         }
         fC = (fC*ANGSTROMS_PER_NM); //-->gaussSigma1
         fC = 1.0/(2 * fC * fC); // Normalizes sigma
         fD = (fD*ANGSTROMS_PER_NM); //-->gaussRepulsive
         pairTable.addPairGaussType2(indexA,indexB,fA,fB,fC,fD);
       } else if (pairFunction == 7) {
         
         // Double well Guassian function
         fA = (fA/JOULES_PER_CALORIE); //-->gaussA
         fB = (fB*ANGSTROMS_PER_NM); //-->gaussMu1
         if(fC == 0 || fE == 0) {
           char buff[100];
           sprintf(buff,"GromacsTopFile.C::Attempting to load zero into gaussSigma.  Please check the pair: %s\n",buf);
           NAMD_die(buff);
         }
         if((fC < 0 || fE < 0)&& !bool_negative_number_warning_flag) {
           iout << iWARN << "Attempting to load a negative standard deviation into the gaussSigma.  Taking the absolute value of the standard deviation.";
           bool_negative_number_warning_flag = true;
         }
         fC = (fC*ANGSTROMS_PER_NM); //-->gaussSigma1
         fC = 1.0/(2 * fC * fC); // Normalizes sigma
         fD = (fD*ANGSTROMS_PER_NM); //-->gaussMu2
         fE = (fE*ANGSTROMS_PER_NM); //-->gaussSigma2
         fE = 1.0/(2 * fE * fE); // Normalizes sigma
         fF = (fE*ANGSTROMS_PER_NM); //-->gaussRepulsive
         pairTable.addPairGaussType2(indexA,indexB,fA,fB,fC,fD,fE,fF);
       } else {
         
         // Generic error statement
         fprintf(stderr,"Unknown pairFunction in GromacsTopFile.C under the PAIRS section: %d\n",pairFunction);
       }
       break;
     case EXCLUSIONS:
       /* Start of JLai modifications August 16th, 2012 */
       if(2 != sscanf(buf," %d %d ",&atomi,&atomj)) {
         fprintf(stderr,"Syntax error in EXCLUSIONS: %s\n",buf);
         exit(1);
       }
       
       // Shift the atom indices to be zero-based
       atomi--;
       atomj--;
 
       /*Load exclusion information into file*/
       exclusions_atom_i.add(atomi);
       exclusions_atom_j.add(atomj);
       numExclusion++;
 
       /* Reading in exclusions information from file and loading */
       break;
       // End of JLai modifications August 16th, 2012
     }
   }
 
   fclose(f);
 }

Member Function Documentation

◆ getAngle()

void GromacsTopFile::getAngle	(	int	num,
		int *	atomi,
		int *	atomj,
		int *	atomk,
		int *	angletype
	)		const

Definition at line 727 of file GromacsTopFile.C.

References GenericAngle::getAtomi(), GenericAngle::getAtomj(), GenericAngle::getAtomk(), GenericAngle::getType(), and ResizeArray< Elem >::size().

                                                 {
   /* 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;
   }
 }

◆ getAngleParams()

void GromacsTopFile::getAngleParams	(	int	num,
		Real *	th0,
		Real *	kth,
		int *	funct
	)		const

Definition at line 1084 of file GromacsTopFile.C.

References AngleTable::getParams().

         {
   angleTable.getParams(num,th0,kth,funct);
 }

◆ getAtom()

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

Definition at line 809 of file GromacsTopFile.C.

References GenericAtom::getAtomName(), GenericAtom::getCharge(), GenericAtom::getMass(), GenericAtom::getResNum(), GenericAtom::getResType(), GenericAtom::getType(), GenericAtom::getTypeNum(), and ResizeArray< Elem >::size().

Referenced by PDB::PDB().

         {
   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;
   }
 }

◆ getAtomParams()

void GromacsTopFile::getAtomParams	(	int	num,
		char *	type
	)		const

inline

Definition at line 544 of file GromacsTopFile.h.

References AtomTable::getType().

Referenced by getVDWParams().

                                                 {
     atomTable.getType(num,type);
   }

◆ getBond()

void GromacsTopFile::getBond	(	int	num,
		int *	atomi,
		int *	atomj,
		int *	bondtype
	)		const

Definition at line 693 of file GromacsTopFile.C.

References GenericBond::getAtomi(), GenericBond::getAtomj(), GenericBond::getType(), and ResizeArray< Elem >::size().

                                                                                {
   /* 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;
   }
 }

◆ getBondParams()

void GromacsTopFile::getBondParams	(	int	num,
		Real *	b0,
		Real *	kB,
		int *	funct
	)		const

Definition at line 1077 of file GromacsTopFile.C.

References BondTable::getParams().

         {
   bondTable.getParams(num,b0,kB,funct);
 }

◆ getDihedral()

void GromacsTopFile::getDihedral	(	int	num,
		int *	atomi,
		int *	atomj,
		int *	atomk,
		int *	atoml,
		int *	type
	)		const

Definition at line 762 of file GromacsTopFile.C.

References GenericDihedral::getAtomi(), GenericDihedral::getAtomj(), GenericDihedral::getAtomk(), GenericDihedral::getAtoml(), GenericDihedral::getType(), and ResizeArray< Elem >::size().

                                                            {
   /* 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;
   }
 }

◆ getDihedralParams()

void GromacsTopFile::getDihedralParams	(	int	num,
		Real *	c,
		int *	mult,
		int *	funct
	)		const

Definition at line 1089 of file GromacsTopFile.C.

References DihedralTable::getParams().

         {
   dihedralTable.getParams(num,c,mult,funct);
 }

◆ getExclusions()

void GromacsTopFile::getExclusions	(	int *	atomi,
		int *	atomj
	)		const

Definition at line 680 of file GromacsTopFile.C.

References exclusions_atom_i, exclusions_atom_j, and ResizeArray< Elem >::size().

                                                                {
     for(int i =0; i < exclusions_atom_i.size(); i++) {
         atomi[i] = exclusions_atom_i[i];
         atomj[i] = exclusions_atom_j[i];
     }
   return;
 }

◆ getNumAngleParams()

int GromacsTopFile::getNumAngleParams ( ) const

inline

Definition at line 522 of file GromacsTopFile.h.

References AngleTable::size().

522 { return angleTable.size(); }

AngleTable::size

int size() const

Definition: GromacsTopFile.h:290

◆ getNumAngles()

int GromacsTopFile::getNumAngles ( ) const

Definition at line 639 of file GromacsTopFile.C.

References ResizeArray< Elem >::size().

                                        {
   int n=0,i;
   for(i=0;i<molInsts.size();i++) {
     n += molInsts[i]->getNum() *
          molInsts[i]->getMol()->getNumAngles();
   }
   return n;
 }

◆ getNumAtomParams()

int GromacsTopFile::getNumAtomParams ( ) const

inline

Definition at line 520 of file GromacsTopFile.h.

References AtomTable::size().

520 { return atomTable.size(); }

AtomTable::size

int size() const

Definition: GromacsTopFile.h:227

◆ getNumAtoms()

int GromacsTopFile::getNumAtoms ( ) const

Definition at line 796 of file GromacsTopFile.C.

References ResizeArray< Elem >::size().

Referenced by PDB::PDB().

                                       {
   int n=0;
   int i;
   for(i=0;i<molInsts.size();i++) {
     n += molInsts[i]->getNum() *
          molInsts[i]->getMol()->getNumAtoms();
   }
   return n;
 }

◆ getNumBondParams()

int GromacsTopFile::getNumBondParams ( ) const

inline

Definition at line 521 of file GromacsTopFile.h.

References BondTable::size().

521 { return bondTable.size(); }

BondTable::size

int size() const

Definition: GromacsTopFile.h:247

◆ getNumBonds()

int GromacsTopFile::getNumBonds ( ) const

Definition at line 629 of file GromacsTopFile.C.

References ResizeArray< Elem >::size().

                                       {
   int n=0,i;
   for(i=0;i<molInsts.size();i++) {
     n += molInsts[i]->getNum() *
          molInsts[i]->getMol()->getNumBonds();
   }
   return n;
 }

◆ getNumDihedralParams()

int GromacsTopFile::getNumDihedralParams ( ) const

inline

Definition at line 523 of file GromacsTopFile.h.

References DihedralTable::size().

523 { return dihedralTable.size(); }

DihedralTable::size

int size() const

Definition: GromacsTopFile.h:349

◆ getNumDihedrals()

int GromacsTopFile::getNumDihedrals ( ) const

Definition at line 649 of file GromacsTopFile.C.

References ResizeArray< Elem >::size().

                                           {
   int n=0,i;
   for(i=0;i<molInsts.size();i++) {
     n += molInsts[i]->getNum() *
          molInsts[i]->getMol()->getNumDihedrals();
   }
   return n;
 }

◆ getNumExclusions()

int GromacsTopFile::getNumExclusions ( ) const

Definition at line 675 of file GromacsTopFile.C.

References numExclusion.

                                            {
   return numExclusion;
 }

◆ getNumGaussPair()

int GromacsTopFile::getNumGaussPair ( ) const

Definition at line 670 of file GromacsTopFile.C.

References numGaussPair.

                                           {
   return numGaussPair;
 }

◆ getNumLJPair()

int GromacsTopFile::getNumLJPair ( ) const

Definition at line 666 of file GromacsTopFile.C.

References numLJPair.

                                        {
   return numLJPair;
 }

◆ getNumPair()

int GromacsTopFile::getNumPair ( ) const

Definition at line 660 of file GromacsTopFile.C.

References numGaussPair, and numLJPair.

                                      {
   int numPair = 0;
   numPair = numLJPair + numGaussPair;
   return numPair;
 }

◆ getPairGaussArrays2()

void GromacsTopFile::getPairGaussArrays2	(	int *	indexA,
		int *	indexB,
		Real *	gaussA,
		Real *	gaussMu1,
		Real *	gaussSigma1,
		Real *	gaussMu2,
		Real *	gaussSigma2,
		Real *	gaussRepulsive
	)

Definition at line 1098 of file GromacsTopFile.C.

References PairTable::getPairGaussArrays2().

                                                           {
   pairTable.getPairGaussArrays2(indexA, indexB, gaussA, gaussMu1, gaussSigma1, gaussMu2, gaussSigma2, gaussRepulsive);
 }

◆ getPairLJArrays2()

void GromacsTopFile::getPairLJArrays2	(	int *	indexA,
		int *	indexB,
		Real *	pairC6,
		Real *	pairC12
	)

Definition at line 1094 of file GromacsTopFile.C.

References PairTable::getPairLJArrays2().

                                                                                            {
   pairTable.getPairLJArrays2(indexA, indexB, pairC6, pairC12);
 }

◆ getSystemName()

char* GromacsTopFile::getSystemName ( ) const

inline

Definition at line 511 of file GromacsTopFile.h.

511 { return systemName; }

◆ getVDWParams()

void GromacsTopFile::getVDWParams	(	int	typea,
		int	typeb,
		Real *	c6,
		Real *	c12,
		Real *	c6pair,
		Real *	c7
	)		const

Definition at line 1417 of file GromacsTopFile.C.

References getAtomParams(), AtomTable::getParams(), and VDWTable::getParams().

                                                                           {
   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) {
     // QQ151069
     /*if ( !genPairs && numa != numb ) {
       iout << iWARN << "VDW table using combining rule for types "
       << typea << " " << typeb << "\n" << endi;
       }*/
 
     /* 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));
 
     /*  this is the only reasonable option  */
     *c6pair  = *c6  * fudgeLJ;
     *c12pair = *c12 * fudgeLJ;
 
   }
 
 }

The documentation for this class was generated from the following files:

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

◆ GromacsTopFile()

Member Function Documentation

◆ getAngle()

◆ getAngleParams()

◆ getAtom()

◆ getAtomParams()

◆ getBond()

◆ getBondParams()

◆ getDihedral()

◆ getDihedralParams()

◆ getExclusions()

◆ getNumAngleParams()

◆ getNumAngles()

◆ getNumAtomParams()

◆ getNumAtoms()

◆ getNumBondParams()

◆ getNumBonds()

◆ getNumDihedralParams()

◆ getNumDihedrals()

◆ getNumExclusions()

◆ getNumGaussPair()

◆ getNumLJPair()

◆ getNumPair()

◆ getPairGaussArrays2()

◆ getPairLJArrays2()

◆ getSystemName()

◆ getVDWParams()