GromacsTopFile Class Reference

#include <GromacsTopFile.h>

List of all members.

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 477 of file GromacsTopFile.h.


Constructor & Destructor Documentation

GromacsTopFile::GromacsTopFile ( char *  filename  ) 

Definition at line 65 of file GromacsTopFile.C.

References ResizeArray< Elem >::add(), PairTable::addPairGaussType2(), PairTable::addPairLJType2(), AtomTable::addType(), DihedralTable::addType(), AngleTable::addType(), BondTable::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, f, DihedralTable::getIndex(), AngleTable::getIndex(), BondTable::getIndex(), AtomTable::getParams(), DihedralTable::getParams(), AngleTable::getParams(), BondTable::getParams(), iout, iWARN(), j, JOULES_PER_CALORIE, LINESIZE, LONGNAMESIZE, MOLECULES, MOLECULETYPE, NAMD_die(), NAMESIZE, NONBOND, numExclusion, PAIRS, PI, ResizeArray< Elem >::size(), SYSTEM, and UNKNOWN.

00065                                              {
00066   fudgeLJ = fudgeQQ = 1.0;
00067   /* open the file */
00068   FILE *f = fopen(filename,"r");
00069   char buf[LINESIZE];
00070   char modename[20];
00071   int mode;
00072   if(f==NULL) {
00073     sprintf(buf,"Error opening file '%s'",filename);
00074     NAMD_die(buf);
00075   }
00076 
00077   /* really bad parser XXX probably just works on the files we
00078      happen to have REWRITE THIS SOON.  It should allow for \- line
00079      continuations, check for errors in the file, etc. */
00080   while(fgets(buf,LINESIZE-1,f)) {
00081     char testchar;
00082     int i,j;
00083 
00084     /* defaults */
00085     int nbfunc, combrule;
00086     char genpairs[20];
00087     
00088     /* atom buffers */
00089     int num, resnum, chargegp, typenum;
00090     char type[NAMESIZE+1], restype[NAMESIZE+1], atomname[NAMESIZE+1];
00091     char particletype[NAMESIZE+1];
00092     float charge, mass, c6, c12, junkf;
00093 
00094     /* moltype buffers */
00095     int nrexcl;
00096     char molname[LONGNAMESIZE+1];
00097 
00098     /* molInst buffers */
00099     int copies;
00100     MolInst *moleculeinstance;
00101 
00102     /* general atomset buffers */
00103     int atomi, atomj, atomk, atoml;
00104     char typea[NAMESIZE+1],typeb[NAMESIZE+1],
00105       typec[NAMESIZE+1],typed[NAMESIZE+1];
00106     const char *tmptypea,*tmptypeb,*tmptypec,*tmptyped;
00107     int funct, index;
00108     float c0,c1;
00109     
00110     /* bonds */
00111     float b0,kB,th0,kth;
00112 
00113     /* dihedrals */
00114     float c[6];
00115     int mult=0;
00116 
00117     /* check for comments */
00118     if(sscanf(buf," %c",&testchar)==1) {
00119       if(testchar == ';') continue;
00120     }
00121     else { /* this is a blank line */
00122       continue;
00123     }
00124 
00125     /* check for a new mode */
00126     if(sscanf(buf," [ %19[^] ] ]",modename)==1) {
00127       /* switch the mode */
00128       if(0==strcmp(modename,"atoms"))              mode = ATOMS;
00129       else if(0==strcmp(modename,"atomtypes"))     mode = ATOMTYPES;
00130       else if(0==strcmp(modename,"moleculetype"))  mode = MOLECULETYPE;
00131       else if(0==strcmp(modename,"molecules"))     mode = MOLECULES;
00132       else if(0==strcmp(modename,"system"))        mode = SYSTEM;
00133       else if(0==strcmp(modename,"bonds"))         mode = BONDS;
00134       else if(0==strcmp(modename,"bondtypes"))     mode = BONDTYPES;
00135       else if(0==strcmp(modename,"angles"))        mode = ANGLES;
00136       else if(0==strcmp(modename,"angletypes"))    mode = ANGLETYPES;
00137       else if(0==strcmp(modename,"dihedrals"))     mode = DIHEDRALS;
00138       else if(0==strcmp(modename,"dihedraltypes")) mode = DIHEDRALTYPES;
00139       else if(0==strcmp(modename,"defaults"))      mode = DEFAULTS;
00140       else if(0==strcmp(modename,"nonbond_params")) mode = NONBOND;
00141       // JLai
00142       else if(0==strcmp(modename,"pairs")) mode = PAIRS;
00143       else if(0==strcmp(modename,"exclusions")) mode = EXCLUSIONS;
00144       else {    
00145         fprintf(stderr,"Warning: unknown mode %s\n",modename);
00146         mode = UNKNOWN;
00147       }
00148 
00149       continue;
00150     }
00151 
00152     /* now do the appropriate thing based on the current mode */
00153     switch(mode) {
00154     case SYSTEM:
00155       systemName = strdup(buf);
00156       break;
00157 
00158     case DEFAULTS:
00159       i = sscanf(buf," %d %d %20s %f %f",
00160                  &nbfunc,&combrule,genpairs,&fudgeLJ,&fudgeQQ);
00161       if(i < 3) { // didn't get enough parameters 
00162         fprintf(stderr,"syntax error in DEFAULTS\n");
00163         exit(1);
00164       }
00165       if(nbfunc != 1) { // I don't know how it works with nbfunc=2
00166         fprintf(stderr,"Non-bonded function != 1 unsupported in DEFAULTS\n");
00167         exit(1);
00168       }
00169       if(combrule != 1) { // same here
00170         fprintf(stderr,"Combination rule != 1 unsupported in DEFAULTS\n");
00171         exit(1);
00172       }
00173       if(0==strcmp(genpairs,"yes")) {
00174         genPairs=1;
00175         if(i!=5) {
00176           fprintf(stderr,"syntax error in DEFAULTS\n");
00177           exit(1);
00178         }
00179         // else fudgeLJ and fudgeQQ got written automatically
00180       }
00181       else genPairs=0;
00182 
00183       break;
00184 
00185     case NONBOND:
00186       if(5 != sscanf(buf," %5s %5s %d %f %f",
00187                      typea, typeb, &funct, &c6, &c12)) {
00188         fprintf(stderr,"Syntax error in NONBOND\n");
00189         exit(1);
00190       }
00191       // convert kJ/mol*nm6 to kcal/mol*A6 and ..12 ..12
00192       c6 =  c6/JOULES_PER_CALORIE*1E6;
00193       c12= c12/JOULES_PER_CALORIE*1E12;
00194       vdwTable.addType(typea,typeb,c6,c12);
00195       break;
00196 
00197     case BONDS:
00198       i = sscanf(buf," %d %d %d %f %f",
00199                  &atomi,&atomj,&funct,&c0,&c1);
00200       atomi--; // shift right away to a zero-indexing
00201       atomj--;
00202       if(i==3) {
00203        tmptypea = genericMols[genericMols.size()-1]->getAtom(atomi)->getType();
00204        tmptypeb = genericMols[genericMols.size()-1]->getAtom(atomj)->getType();
00205         /* find the index and parameters */
00206         index = bondTable.getParams(tmptypea, tmptypeb, funct, &b0, &kB);
00207         if(index==-1) {
00208           fprintf(stderr,"Required bondtype %s--%s (function %d) not found.\n",
00209                   tmptypea,tmptypeb,funct);
00210           exit(1);
00211         }
00212       }
00213       else if(i==5) {
00214         /* first set the values of b0 and kB correctly */
00215         b0 = c0*ANGSTROMS_PER_NM; /* convert nm to A */
00216         if(funct==1) { /* harmonic potential */
00217           /* convert kJ/nm2 to kcal/A2 and use E=kx2 instead of half that. */
00218           kB = c1/JOULES_PER_CALORIE/100/2;
00219         }
00220         else if(funct==2) { /* special fourth-order potential */
00221           /* convert to the normal harmonic constant and kJ/nm2 to kcal/A2 */
00222           kB = 2*c1*c0*c0/JOULES_PER_CALORIE/100;
00223           kB /= 2; /* use the NAMD system where E=kx2 */
00224         }
00225         else {
00226           fprintf(stderr,"I don't know what funct=%d means in BONDS\n",funct);
00227           exit(1);
00228         }
00229         /* look up the index */
00230         index = bondTable.getIndex(b0,kB,funct);
00231       }
00232       else {
00233         fprintf(stderr,"Syntax error in BONDS\n");
00234         exit(1);
00235       }
00236 
00237       genericMols[genericMols.size()-1]->addBond(atomi,atomj,index);
00238       break;
00239       
00240     case BONDTYPES:
00241       if(5 != sscanf(buf," %5s %5s %d %f %f",
00242                      typea,typeb,&funct,&c0,&c1)) {
00243         fprintf(stderr,"Syntax error in BONDTYPES\n");
00244         exit(1);
00245       }
00246 
00247       /* first set the values of b0 and kB correctly */
00248       b0 = c0*10; /* convert nm to A */
00249       if(funct==1) { /* harmonic potential */
00250         /* convert kJ/nm2 to kcal/A2 and use E=kx2 instead of half that. */
00251         kB = c1/JOULES_PER_CALORIE/100/2;
00252       }
00253       else if(funct==2) { /* special fourth-order potential */
00254         /* convert to the normal harmonic constant and kJ/nm2 to kcal/A2 */
00255         kB = 2*c1*c0*c0/JOULES_PER_CALORIE/100;
00256         kB /= 2; /* use the NAMD system where E=kx2 */
00257       }
00258       else {
00259  fprintf(stderr,"I don't know what funct=%d means in BONDTYPES\n",funct);
00260         exit(1);
00261       }
00262 
00263       bondTable.addType(typea,typeb,b0,kB,funct);
00264       break;
00265 
00266     case ANGLES:
00267       i = sscanf(buf," %d %d %d %d %f %f",
00268                      &atomi,&atomj,&atomk,&funct,&c0,&c1);
00269       atomi--; // shift right away to a zero-indexing
00270       atomj--;
00271       atomk--;
00272       if(i == 4) { /* we have to look up the last two parameters */
00273        tmptypea = genericMols[genericMols.size()-1]->getAtom(atomi)->getType();
00274        tmptypeb = genericMols[genericMols.size()-1]->getAtom(atomj)->getType();
00275        tmptypec = genericMols[genericMols.size()-1]->getAtom(atomk)->getType();
00276         /* find the index and parameters */
00277         index = angleTable.getParams(tmptypea, tmptypeb, tmptypec,
00278                                         funct, &th0, &kth);
00279         if(index==-1) {
00280           fprintf(stderr,
00281                   "Required angletype %s--%s--%s (function %d) not found.\n",
00282                   tmptypea,tmptypeb,tmptypec,funct);
00283           exit(1);
00284         }
00285       }
00286       else if(i == 6) {
00287         /* first set the values of th0 and kth correctly */
00288         if(funct == 1) {
00289           th0 = c0; /* both are in degrees */
00290           kth = c1/JOULES_PER_CALORIE/2; /* convert kJ/rad2 to kcal/rad2 and use E=kx2 */
00291         }
00292         else if(funct == 2) {
00293           th0 = c0; /* both are in degrees */
00294           /* convert G96 kJ to kcal/rad2 and use E=kx2 */
00295           kth = sin(th0*PI/180)*sin(th0*PI/180)*c1/JOULES_PER_CALORIE/2;
00296         }
00297         else {
00298           fprintf(stderr,"I don't know what funct=%d means in ANGLES\n",funct);
00299           exit(1);
00300         }
00301         /* add the angle type to our table */
00302         index = angleTable.getIndex(th0,kth,funct);
00303       }
00304       else {
00305         fprintf(stderr,"Syntax error (%d args) in ANGLES: %s\n",i,buf);
00306         exit(1);
00307       }
00308 
00309       /* add the angle to our table */
00310       genericMols[genericMols.size()-1]->addAngle(atomi,atomj,atomk,index);
00311       break;
00312 
00313     case ANGLETYPES:
00314       if(6 != sscanf(buf," %5s %5s %5s %d %f %f",
00315                      typea,typeb,typec,&funct,&c0,&c1)) {
00316         fprintf(stderr,"Syntax error in ANGLETYPES\n");
00317         exit(1);
00318       }
00319       /* first set the values of th0 and kth correctly */
00320       if(funct == 1) {
00321         th0 = c0; /* both are in degrees */
00322         kth = c1/JOULES_PER_CALORIE/2; /* convert kJ/rad2 to kcal/rad2 and use E=kx2 */
00323       }
00324       else if(funct == 2) {
00325         th0 = c0; /* both are in degrees */
00326         /* convert G96 kJ to kcal/rad2 and use E=kx2 */
00327         kth = sin(th0*PI/180)*sin(th0*PI/180)*c1/JOULES_PER_CALORIE/2;
00328       }
00329       else {
00330         fprintf(stderr,"I don't know what funct=%d means in ANGLETYPES\n",
00331                 funct);
00332         exit(1);
00333       }
00334       angleTable.addType(typea,typeb,typec,th0,kth,funct);
00335       break;
00336 
00337     case DIHEDRALS:
00338       i = sscanf(buf," %d %d %d %d %d %f %f %f %f %f %f",
00339                  &atomi,&atomj,&atomk,&atoml,&funct,
00340                  &c[0],&c[1],&c[2],&c[3],&c[4],&c[5]);
00341       atomi--; // shift right away to a zero-indexing
00342       atomj--;
00343       atomk--;
00344       atoml--;
00345       if(i==5) { /* we have to look up the parameters */
00346        tmptypea = genericMols[genericMols.size()-1]->getAtom(atomi)->getType();
00347        tmptypeb = genericMols[genericMols.size()-1]->getAtom(atomj)->getType();
00348        tmptypec = genericMols[genericMols.size()-1]->getAtom(atomk)->getType();
00349        tmptyped = genericMols[genericMols.size()-1]->getAtom(atoml)->getType();
00350         /* find the index and parameters */
00351         index = dihedralTable.getParams(tmptypea, tmptypeb, tmptypec,
00352                                         tmptyped, funct, c, &mult);
00353         if(index==-1) {
00354           fprintf(stderr,
00355              "Required dihedraltype %s--%s--%s--%s (function %d) not found.\n",
00356                   tmptypea,tmptypeb,tmptypec,tmptyped,funct);
00357           exit(1);
00358         }
00359       }
00360       else if(i==7 || i==8 || i==11) { /* the parameters are given */
00361         if(funct==1 || funct==2) { /* we should have two parameters */
00362           if(i!=7+(funct==1)) {    /* plus a multiplicity for funct==1 */
00363             fprintf(stderr,"Must have 7 args for funct=1,2\n");
00364             exit(1);
00365           }
00366           c[0] = c[0]; /* both in deg */
00367           if(i==7) {
00368               c[1] = c[1]/(2*JOULES_PER_CALORIE); /* convert kJ to kcal and still use E=kx2*/
00369           } else if (i==8 || i==11) {
00370               c[1] = c[1]/(1*JOULES_PER_CALORIE); /* convert kJ to kcal and still use E=kx2*/
00371           } 
00372           //c[1] = c[1]/(2*JOULES_PER_CALORIE); /* convert kJ to kcal and still use E=kx2*/
00373           /* for funct==1 these are both divided by rad^2 */
00374           if(funct==1) {
00375             mult=(int)(c[2]+0.5); /* round to nearest integer :p */
00376           }
00377         }
00378         else if(funct==3) {
00379           if(i!=11){
00380             fprintf(stderr,"Must have 11 args for funct=3\n");
00381             exit(1);
00382           }
00383 
00384           for(j=0;j<6;j++) {
00385             c[j] = c[j]/JOULES_PER_CALORIE; /* convert kJ to kcal and E=Cn cos^n */
00386           }
00387         }
00388         else {
00389           fprintf(stderr,
00390                   "I don't know what funct=%d means in DIHEDRALS\n",funct);
00391           exit(1);
00392         }
00393         index = dihedralTable.getIndex(c,mult,funct);
00394       }
00395       else {
00396         fprintf(stderr,"Syntax error (%d args) in DIHEDRALS: %s\n",i,buf);
00397         exit(1);
00398       }
00399 
00400       /* add the dihedrals to our table */
00401       genericMols[genericMols.size()-1]->addDihedral(atomi,atomj,atomk,atoml,
00402                                                      index);
00403       break;
00404     case DIHEDRALTYPES:
00405       i = sscanf(buf," %5s %5s %d %f %f %f %f %f %f",
00406                  typea,typeb,&funct,&c[0],&c[1],&c[2],&c[3],&c[4],&c[5]);
00407       if(funct == 1 || funct == 2) {
00408         if(i!=5+(funct==1)) { /* 6 for f=2, 5 for f=1 */
00409           fprintf(stderr,"Syntax error in DIHEDRALTYPES: %s\n",buf);
00410           exit(1);
00411         }
00412         c[0] = c[0]; /* both in deg */
00413         c[1] = c[1]/JOULES_PER_CALORIE; /* convert kJ to kcal and still use E=kx2*/
00414         /* for funct==1 these are both divided by rad^2 */
00415         if(funct==1) {
00416           mult=(int)(c[2]+0.5); /* round to nearest integer :p */
00417         }
00418       }
00419       else if(funct == 3) {
00420         if(i!=9) {
00421           fprintf(stderr,"Syntax error in DIHEDRALTYPES\n");
00422           exit(1);
00423         }
00424         for(j=0;j<6;j++) {
00425           c[j] = c[j]/JOULES_PER_CALORIE; /* convert kJ to kcal and E=Cn cos^n */
00426         }
00427       }
00428       else {
00429         fprintf(stderr,"I don't know what funct=%d means in DIHEDRALTYPES\n",
00430                 funct);
00431         exit(1);
00432       }
00433       dihedralTable.addType(typea,typeb,c,mult,funct);
00434       break;
00435       
00436     case ATOMS:
00437       i = sscanf(buf," %d %5s %d %5s %5s %d %f %f",
00438                    &num, type, &resnum, restype,
00439                  atomname, &chargegp, &charge, &mass);
00440       if(i==7) { /* XXX temporary - I should be able to get more
00441                     params */
00442         typenum = atomTable.getParams(type,&mass,&junkf,&junkf,&junkf);
00443         i=8;
00444       }
00445       else {
00446         if(i!=8) {
00447           fprintf(stderr,"Syntax error in ATOMS\n");
00448           exit(1);
00449         }
00450         // just get the type number
00451         typenum = atomTable.getParams(type,&junkf,&junkf,&junkf,&junkf);
00452       }
00453       genericMols[genericMols.size()-1]->addAtom(type,typenum,resnum,
00454                                                  restype,atomname,charge,mass);
00455       break;
00456 
00457     case ATOMTYPES:
00458       if(6 != sscanf(buf," %5s %f %f %5s %f %f",type,&mass,&charge,
00459                      particletype,&c6,&c12)) {
00460         fprintf(stderr,"Syntax error in ATOMTYPES: %s\n",buf);
00461         exit(1);
00462       }
00463       /* conversions:
00464          c6  - kJ/mol nm6  -> kcal/mol A6
00465          c12 - kJ/mol nm12 -> kcal/mol A12 */
00466       atomTable.addType(type,mass,charge,
00467                         c6/(JOULES_PER_CALORIE)*1E6,
00468                         c12/(JOULES_PER_CALORIE)*1E12);
00469       break;
00470 
00471     case MOLECULETYPE:
00472       if(2!=sscanf(buf," %20s %d",molname,&nrexcl)) {
00473         fprintf(stderr,"Syntax error in MOLECULETYPE: %s\n",buf);
00474         exit(1);
00475       }
00476 
00477       /* add another generic molecule holder */
00478       genericMols.add(new GenericMol(molname));
00479       break;
00480 
00481     case MOLECULES:
00482       if(2!=sscanf(buf," %20s %d",molname,&copies)) {
00483         fprintf(stderr,"Syntax error in MOLECULES: %s\n",buf);
00484         exit(1);
00485       }
00486       
00487       /* search for the specified molecule and make a molInst of it*/
00488       moleculeinstance = NULL;
00489       for(i=0;i<genericMols.size();i++) {
00490         if(0==strcmp(molname,genericMols[i]->getName())) {
00491           moleculeinstance = new MolInst(genericMols[i],copies);
00492           break;
00493         }
00494       }
00495 
00496       if(moleculeinstance==NULL) {
00497         fprintf(stderr,"Molecule %s undefined.\n",molname);
00498         exit(1);
00499       }
00500       
00501       /* put it at the end of the list */
00502       molInsts.add(moleculeinstance);
00503 
00504       break;
00505     case PAIRS:
00506       int indexA;
00507       int indexB;
00508       int pairFunction;
00509       Real fA;
00510       Real fB;
00511       Real fC;
00512       Real fD;
00513       Real fE;
00514       Real fF;
00515       
00516       int countVariables;
00517       countVariables = sscanf(buf," %d %d %d %f %f %f %f %f %f",&indexA,&indexB,&pairFunction,&fA,&fB,&fC,&fD,&fE,&fF);
00518 
00519       if ((countVariables >= 4 && countVariables >= 10)) {
00520         fprintf(stderr,"Syntax error in PAIRS: %s\n",buf);
00521         exit(1);
00522       }
00523       
00524       // Shift the atom indices to be zero-based
00525       indexA--;
00526       indexB--;
00527 
00528       // Make sure that the indexA is less than indexB
00529       /*if ( indexA > indexB ) {
00530           int tmpIndex = indexA;
00531           indexB = indexA;
00532           indexA = tmpIndex;
00533           }*/
00534 
00535       if (pairFunction == 1) {
00536         
00537         // LJ code
00538         fA = (fA/JOULES_PER_CALORIE)*1E6;
00539         fB= (fB/JOULES_PER_CALORIE)*1E12;
00540         pairTable.addPairLJType2(indexA,indexB,fA,fB);
00541       } else if (pairFunction == 5){
00542         
00543         // Bare Gaussian potential
00544         fA = (fA/JOULES_PER_CALORIE); //-->gaussA
00545         fB = (fB*ANGSTROMS_PER_NM); //-->gaussMu1
00546         if(fC == 0) {
00547           char buff[100];
00548           sprintf(buff,"GromacsTopFile.C::Attempting to load zero into gaussSigma.  Please check the pair: %s\n",buf);
00549           NAMD_die(buff);
00550         }
00551         if(fC < 0 && !bool_negative_number_warning_flag) {
00552           iout << iWARN << "Attempting to load a negative standard deviation into the gaussSigma.  Taking the absolute value of the standard deviation.";
00553           bool_negative_number_warning_flag = true;
00554         }
00555         fC = (fC*ANGSTROMS_PER_NM); //-->gaussSigma1
00556         fC = 1.0/(2 * fC * fC); // Normalizes sigma
00557         pairTable.addPairGaussType2(indexA,indexB,fA,fB,fC);
00558       } else if (pairFunction == 6) {
00559         
00560         // Combined Gaussian + repulsive r^-12 potential
00561         fA = (fA/JOULES_PER_CALORIE); //-->gaussA
00562         fB = (fB*ANGSTROMS_PER_NM); //-->gaussMu1
00563         if(fC == 0) {
00564           char buff[100];
00565           sprintf(buff,"GromacsTopFile.C::Attempting to load zero into gaussSigma.  Please check the pair: %s\n",buf);
00566           NAMD_die(buff);
00567         }
00568         if(fC < 0 && !bool_negative_number_warning_flag) {
00569           iout << iWARN << "Attempting to load a negative standard deviation into the gaussSigma.  Taking the absolute value of the standard deviation.";
00570           bool_negative_number_warning_flag = true;
00571         }
00572         fC = (fC*ANGSTROMS_PER_NM); //-->gaussSigma1
00573         fC = 1.0/(2 * fC * fC); // Normalizes sigma
00574         fD = (fD*ANGSTROMS_PER_NM); //-->gaussRepulsive
00575         pairTable.addPairGaussType2(indexA,indexB,fA,fB,fC,fD);
00576       } else if (pairFunction == 7) {
00577         
00578         // Double well Guassian function
00579         fA = (fA/JOULES_PER_CALORIE); //-->gaussA
00580         fB = (fB*ANGSTROMS_PER_NM); //-->gaussMu1
00581         if(fC == 0 || fE == 0) {
00582           char buff[100];
00583           sprintf(buff,"GromacsTopFile.C::Attempting to load zero into gaussSigma.  Please check the pair: %s\n",buf);
00584           NAMD_die(buff);
00585         }
00586         if((fC < 0 || fE < 0)&& !bool_negative_number_warning_flag) {
00587           iout << iWARN << "Attempting to load a negative standard deviation into the gaussSigma.  Taking the absolute value of the standard deviation.";
00588           bool_negative_number_warning_flag = true;
00589         }
00590         fC = (fC*ANGSTROMS_PER_NM); //-->gaussSigma1
00591         fC = 1.0/(2 * fC * fC); // Normalizes sigma
00592         fD = (fD*ANGSTROMS_PER_NM); //-->gaussMu2
00593         fE = (fE*ANGSTROMS_PER_NM); //-->gaussSigma2
00594         fE = 1.0/(2 * fE * fE); // Normalizes sigma
00595         fF = (fE*ANGSTROMS_PER_NM); //-->gaussRepulsive
00596         pairTable.addPairGaussType2(indexA,indexB,fA,fB,fC,fD,fE,fF);
00597       } else {
00598         
00599         // Generic error statement
00600         fprintf(stderr,"Unknown pairFunction in GromacsTopFile.C under the PAIRS section: %d\n",pairFunction);
00601       }
00602       break;
00603     case EXCLUSIONS:
00604       /* Start of JLai modifications August 16th, 2012 */
00605       if(2 != sscanf(buf," %d %d ",&atomi,&atomj)) {
00606         fprintf(stderr,"Syntax error in EXCLUSIONS: %s\n",buf);
00607         exit(1);
00608       }
00609       
00610       // Shift the atom indices to be zero-based
00611       atomi--;
00612       atomj--;
00613 
00614       /*Load exclusion information into file*/
00615       exclusions_atom_i.add(atomi);
00616       exclusions_atom_j.add(atomj);
00617       numExclusion++;
00618 
00619       /* Reading in exclusions information from file and loading */
00620       break;
00621       // End of JLai modifications August 16th, 2012
00622     }
00623   }
00624 
00625   fclose(f);
00626 }


Member Function Documentation

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().

00728                                                 {
00729   /* figure out what instance this angle is in */
00730   int nangles=0;
00731   int natoms=0;
00732   int i;
00733   for(i=0;i<molInsts.size();i++) {
00734     int molangles = molInsts[i]->getMol()->getNumAngles();
00735     int molatoms = molInsts[i]->getMol()->getNumAtoms();
00736     int newangles = molInsts[i]->getNumAngles();
00737     int newatoms = molInsts[i]->getNumAtoms();
00738     
00739     if(nangles+newangles-1 >= n) {
00740       /* the angle is in this MolInst */
00741       int localnum = (n-nangles) % molangles; /* number within this inst */
00742       int instnum = (n-nangles) / molangles;  /* number of instances before */
00743       int addatoms = natoms+instnum*molatoms; /* extra atms to add to atmi */
00744 
00745       const GenericAngle *a = molInsts[i]->getMol()->getAngle(localnum);
00746 
00747       *atomi = a->getAtomi() + addatoms;
00748       *atomj = a->getAtomj() + addatoms;
00749       *atomk = a->getAtomk() + addatoms;
00750       *angletype = a->getType();
00751       break;
00752     }
00753 
00754     nangles += newangles;
00755     natoms += newatoms;
00756   }
00757 }

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

Definition at line 1084 of file GromacsTopFile.C.

References AngleTable::getParams().

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

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(), getNumAtoms(), GenericAtom::getResNum(), GenericAtom::getResType(), GenericAtom::getType(), GenericAtom::getTypeNum(), and ResizeArray< Elem >::size().

Referenced by PDB::PDB().

00813         {
00814   int natoms=0,n=num; // zero-indexed array
00815   int resnum=0;
00816   /* figure out what instance the atom is in, and what residue number
00817      it has */
00818   int i;
00819 
00820   for(i=0;i<molInsts.size();i++) {
00821     int numnew = molInsts[i]->getNumAtoms(); /* atoms in this MolInst */
00822     int resmol = molInsts[i]->getMol()->getNumRes(); /* # res/mol */
00823     int newres = molInsts[i]->getNumRes(); /* residues in this MolInst */
00824 
00825     if(natoms+numnew-1 >= n) {
00826       /* the atom is in this molInst */
00827       int localnum = (n-natoms) % molInsts[i]->getMol()->getNumAtoms();
00828       int instnum  = (n-natoms) / molInsts[i]->getMol()->getNumAtoms();
00829       
00830       // getAtom is zero-indexed
00831       const GenericAtom *a = molInsts[i]->getMol()->getAtom(localnum);
00832       int residue = resnum + resmol*instnum + a->getResNum();
00833 
00834       *residue_number = residue;
00835       strncpy(residue_name,a->getResType(),11);
00836       strncpy(atom_name,a->getAtomName(),11);
00837       strncpy(atom_type,a->getType(),11);
00838       *charge=a->getCharge();
00839       *mass=a->getMass();
00840       *atom_typenum=a->getTypeNum();
00841       break;
00842     }
00843 
00844     /* go to the next molInst */
00845     natoms += numnew;
00846     resnum += newres;
00847   }
00848 }

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

Definition at line 542 of file GromacsTopFile.h.

References AtomTable::getType().

Referenced by getVDWParams().

00542                                                 {
00543     atomTable.getType(num,type);
00544   }

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().

00693                                                                                {
00694   /* figure out what instance this bond is in */
00695   int nbonds=0;
00696   int natoms=0;
00697   int i;
00698   for(i=0;i<molInsts.size();i++) {
00699     int molbonds = molInsts[i]->getMol()->getNumBonds();
00700     int molatoms = molInsts[i]->getMol()->getNumAtoms();
00701     int newbonds = molInsts[i]->getNumBonds();
00702     int newatoms = molInsts[i]->getNumAtoms();
00703     
00704     if(nbonds+newbonds-1 >= n) {
00705       /* the bond is in this MolInst */
00706       int localnum = (n-nbonds) % molbonds;   /* number within this inst */
00707       int instnum = (n-nbonds) / molbonds;    /* number of instances before */
00708       int addatoms = natoms+instnum*molatoms; /* extra atoms to add to atomi */
00709 
00710       const GenericBond *b = molInsts[i]->getMol()->getBond(localnum);
00711 
00712       *atomi = b->getAtomi() + addatoms;
00713       *atomj = b->getAtomj() + addatoms;
00714       *bondtype = b->getType();
00715       break;
00716     }
00717 
00718     nbonds += newbonds;
00719     natoms += newatoms;
00720   }
00721 }

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

Definition at line 1077 of file GromacsTopFile.C.

References BondTable::getParams().

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

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().

00763                                                            {
00764   /* figure out what instance this angle is in */
00765   int ndihedrals=0;
00766   int natoms=0;
00767   int i;
00768   for(i=0;i<molInsts.size();i++) {
00769     int moldihedrals = molInsts[i]->getMol()->getNumDihedrals();
00770     int molatoms = molInsts[i]->getMol()->getNumAtoms();
00771     int newdihedrals = molInsts[i]->getNumDihedrals();
00772     int newatoms = molInsts[i]->getNumAtoms();
00773     
00774     if(ndihedrals+newdihedrals-1 >= n) {
00775       /* the dihedral is in this MolInst */
00776       int localnum = (n-ndihedrals) % moldihedrals; /* number in this inst */
00777       int instnum = (n-ndihedrals) / moldihedrals; /* number of insts before */
00778       int addatoms = natoms+instnum*molatoms; /*extra atms to add to atmi*/
00779 
00780       const GenericDihedral *a = molInsts[i]->getMol()->getDihedral(localnum);
00781 
00782       *atomi = a->getAtomi() + addatoms;
00783       *atomj = a->getAtomj() + addatoms;
00784       *atomk = a->getAtomk() + addatoms;
00785       *atoml = a->getAtoml() + addatoms;
00786       *type = a->getType();
00787       break;
00788     }
00789 
00790     ndihedrals += newdihedrals;
00791     natoms += newatoms;
00792   }
00793 }

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

Definition at line 1089 of file GromacsTopFile.C.

References DihedralTable::getParams().

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

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

Definition at line 680 of file GromacsTopFile.C.

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

00680                                                                {
00681     for(int i =0; i < exclusions_atom_i.size(); i++) {
00682         atomi[i] = exclusions_atom_i[i];
00683         atomj[i] = exclusions_atom_j[i];
00684     }
00685   return;
00686 }

int GromacsTopFile::getNumAngleParams (  )  const [inline]

Definition at line 520 of file GromacsTopFile.h.

References AngleTable::size().

00520 { return angleTable.size(); }

int GromacsTopFile::getNumAngles (  )  const

Definition at line 639 of file GromacsTopFile.C.

References ResizeArray< Elem >::size().

00639                                        {
00640   int n=0,i;
00641   for(i=0;i<molInsts.size();i++) {
00642     n += molInsts[i]->getNum() *
00643          molInsts[i]->getMol()->getNumAngles();
00644   }
00645   return n;
00646 }

int GromacsTopFile::getNumAtomParams (  )  const [inline]

Definition at line 518 of file GromacsTopFile.h.

References AtomTable::size().

00518 { return atomTable.size(); }

int GromacsTopFile::getNumAtoms (  )  const

Definition at line 796 of file GromacsTopFile.C.

References ResizeArray< Elem >::size().

Referenced by getAtom(), and PDB::PDB().

00796                                       {
00797   int n=0;
00798   int i;
00799   for(i=0;i<molInsts.size();i++) {
00800     n += molInsts[i]->getNum() *
00801          molInsts[i]->getMol()->getNumAtoms();
00802   }
00803   return n;
00804 }

int GromacsTopFile::getNumBondParams (  )  const [inline]

Definition at line 519 of file GromacsTopFile.h.

References BondTable::size().

00519 { return bondTable.size(); }

int GromacsTopFile::getNumBonds (  )  const

Definition at line 629 of file GromacsTopFile.C.

References ResizeArray< Elem >::size().

00629                                       {
00630   int n=0,i;
00631   for(i=0;i<molInsts.size();i++) {
00632     n += molInsts[i]->getNum() *
00633          molInsts[i]->getMol()->getNumBonds();
00634   }
00635   return n;
00636 }

int GromacsTopFile::getNumDihedralParams (  )  const [inline]

Definition at line 521 of file GromacsTopFile.h.

References DihedralTable::size().

00521 { return dihedralTable.size(); }

int GromacsTopFile::getNumDihedrals (  )  const

Definition at line 649 of file GromacsTopFile.C.

References ResizeArray< Elem >::size().

00649                                           {
00650   int n=0,i;
00651   for(i=0;i<molInsts.size();i++) {
00652     n += molInsts[i]->getNum() *
00653          molInsts[i]->getMol()->getNumDihedrals();
00654   }
00655   return n;
00656 }

int GromacsTopFile::getNumExclusions (  )  const

Definition at line 675 of file GromacsTopFile.C.

References numExclusion.

00675                                            {
00676   return numExclusion;
00677 }

int GromacsTopFile::getNumGaussPair (  )  const

Definition at line 670 of file GromacsTopFile.C.

References numGaussPair.

00670                                           {
00671   return numGaussPair;
00672 }

int GromacsTopFile::getNumLJPair (  )  const

Definition at line 666 of file GromacsTopFile.C.

References numLJPair.

00666                                        {
00667   return numLJPair;
00668 }

int GromacsTopFile::getNumPair (  )  const

Definition at line 660 of file GromacsTopFile.C.

References numGaussPair, and numLJPair.

00660                                      {
00661   int numPair = 0;
00662   numPair = numLJPair + numGaussPair;
00663   return numPair;
00664 }

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().

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

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

Definition at line 1094 of file GromacsTopFile.C.

References PairTable::getPairLJArrays2().

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

char* GromacsTopFile::getSystemName (  )  const [inline]

Definition at line 509 of file GromacsTopFile.h.

00509 { return systemName; }

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().

01418                                                                           {
01419   int i,ret;
01420   Real c6a,c6b,c12a,c12b, mjunk, qjunk;
01421   char typea[11]="",typeb[11]="";
01422 
01423   /* get the string names corresponding to numa and numb */
01424   getAtomParams(numa,typea);
01425   getAtomParams(numb,typeb);
01426   if(typea[0]==0) { /* found a bug in my use of strncpy here once */
01427     fprintf(stderr,"Failed to get name of atom %d\n",numa);
01428     exit(1);
01429   }
01430   if(typeb[0]==0) {
01431     fprintf(stderr,"Failed to get name of atom %d\n",numb);
01432     exit(1);
01433   }
01434 
01435   /* first try - use the VDW table */
01436   i = vdwTable.getParams(typea,typeb,c6,c12,c6pair,c12pair);
01437 
01438   if(i==-1) {
01439     // QQ151069
01440     /*if ( !genPairs && numa != numb ) {
01441       iout << iWARN << "VDW table using combining rule for types "
01442       << typea << " " << typeb << "\n" << endi;
01443       }*/
01444 
01445     /* fallback - use the individual atom's parameters */
01446     ret = atomTable.getParams(typea, &mjunk, &qjunk, &c6a, &c12a);
01447     if(ret==-1) {
01448       fprintf(stderr,"Couldn't find atom type %s\n",typea);
01449       exit(1);
01450     }
01451     ret = atomTable.getParams(typeb, &mjunk, &qjunk, &c6b, &c12b);
01452     if(ret==-1) {
01453       fprintf(stderr,"Couldn't find atom type %s\n",typeb);
01454       exit(1);
01455     }
01456     
01457     *c6  = (float)(sqrt(c6a * c6b));
01458     *c12 = (float)(sqrt(c12a*c12b));
01459 
01460     /*  this is the only reasonable option  */
01461     *c6pair  = *c6  * fudgeLJ;
01462     *c12pair = *c12 * fudgeLJ;
01463 
01464   }
01465 
01466 }


The documentation for this class was generated from the following files:
Generated on Thu Nov 23 01:17:18 2017 for NAMD by  doxygen 1.4.7