GromacsTopFile.C

Go to the documentation of this file.
00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <string.h>
00004 #include <math.h>
00005 #include <algorithm>
00006 #include "common.h"
00007 #include "ResizeArray.h"
00008 #include "InfoStream.h"
00009 #include "GromacsTopFile.h"
00010 #include "InfoStream.h"
00011 
00012 #define JOULES_PER_CALORIE 4.184
00013 #define ANGSTROMS_PER_NM 10
00014 
00015 /* A GromacsTopFile represents the information stored in a GROMACS
00016    topolgy file.  This is an immutable type. */
00017 
00018 /* XXX warning: this code contains a few algorithms which run in a
00019    time of O(n^3) or so in the size of the topology file, since I
00020    haven't bothered to do any sorting.  I don't think it matters much,
00021    since it still manages to run on the biggest files in less than a
00022    second, and nothing (memory or time) depends on the total number of
00023    atoms in the simulation - all that matters is the number that have
00024    to be individually defined. */
00025 
00026 /* GromacsTopFile initializes this to represent the data stored in the
00027    file <filename>, or exits on error. */
00028 #define LINESIZE 200
00029 
00030 /* modes */
00031 #define UNKNOWN 0
00032 #define ATOMS 1
00033 #define MOLECULETYPE 2
00034 #define MOLECULES 3
00035 #define SYSTEM 4
00036 #define BONDS 5
00037 #define BONDTYPES 6
00038 #define ANGLES 7
00039 #define ATOMTYPES 8
00040 #define ANGLETYPES 9
00041 #define DIHEDRALS 10
00042 #define DIHEDRALTYPES 11
00043 #define DEFAULTS 12
00044 #define NONBOND 13
00045 #define PAIRS 14
00046 #define EXCLUSIONS 15
00047 #ifndef CODE_REDUNDANT
00048 #define CODE_REDUNDANT 0
00049 #endif
00050 
00051 #define MIN_DEBUG_LEVEL 3
00052 #include "Debug.h"
00053 
00054 /* Initialize variables for exclusion calculation */
00055 /* JLai */
00056 ResizeArray<int> exclusions_atom_i;
00057 ResizeArray<int> exclusions_atom_j;
00058 int numExclusion = 0;
00059 int numLJPair = 0;
00060 int numGaussPair = 0;
00061 bool bool_negative_number_warning_flag = false;
00062 
00063 
00064 
00065 GromacsTopFile::GromacsTopFile(char *filename) {
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 }
00627 
00628 /* returns the number of bonds in the file */
00629 int GromacsTopFile::getNumBonds() const {
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 }
00637 
00638 /* returns the number of angles in the file */
00639 int GromacsTopFile::getNumAngles() const {
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 }
00647 
00648 /* returns the number of dihedral angles in the file */
00649 int GromacsTopFile::getNumDihedrals() const {
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 }
00657 
00658 /* JLai -- August 16th, 2012 modifications*/
00659 /* returns the number of pair interactions in the file */
00660 int GromacsTopFile::getNumPair() const {
00661   int numPair = 0;
00662   numPair = numLJPair + numGaussPair;
00663   return numPair;
00664 }
00665 
00666 int GromacsTopFile::getNumLJPair() const {
00667   return numLJPair;
00668 }
00669 
00670 int GromacsTopFile::getNumGaussPair() const {
00671   return numGaussPair;
00672 }
00673 
00674 /* return the number of exclusions in the file */
00675 int GromacsTopFile::getNumExclusions() const {
00676   return numExclusion;
00677 }
00678 
00679 /* return the list of exclusions from the file */
00680 void GromacsTopFile::getExclusions(int* atomi, int* atomj) const {
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 }
00687 /* End of JLai modifications */
00688 
00689 /* getBond puts the information about bond number <num> into the
00690    spaces pointed to by the other arguments.  Bond number 0 is the
00691    first bond in the list.
00692    For atomi and atomj, 1 is the first atom in the list. */
00693 void GromacsTopFile::getBond(int n, int *atomi, int *atomj, int *bondtype) const {
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 }
00722 
00723 /* getAngle puts the information about angle number <num> into the
00724    spaces pointed to by the other arguments.  Angle number 0 is the
00725    first angle in the list.
00726 */
00727 void GromacsTopFile::getAngle(int n, int *atomi, int *atomj, int *atomk,
00728                           int *angletype) const {
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 }
00758 
00759 /* getDihedral puts the information about dihedral number <num> into
00760    the spaces pointed to by the other arguments.  Dihedral number 0
00761    is the first angle in the list. */
00762 void GromacsTopFile::getDihedral(int n, int *atomi, int *atomj, int *atomk,
00763                               int *atoml, int *type) const {
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 }
00794 
00795 /* getNumAtoms returns the number of atoms stored in the file. */
00796 int GromacsTopFile::getNumAtoms() const {
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 }
00805 
00806 /* getAtom puts the information about atom number <n> into the
00807    spaces pointed to by the other arguments.  The string buffers
00808    must be at least 11 characters long. */
00809 void GromacsTopFile::getAtom(int num,
00810                int *residue_number, char *residue_name,
00811                char *atom_name, char *atom_type, int *atom_typenum,
00812                Real *charge, Real *mass) 
00813   const {
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 }
00849 
00850 GenericAtom::GenericAtom(const char *theType, int theTypeNum, int theResNum,
00851                          const char *theResType,const char *theAtomName,
00852                          Real theCharge, Real theMass) {
00853   strncpy(type,theType,NAMESIZE+1);
00854   typeNum = theTypeNum;
00855   resNum = theResNum;
00856   strncpy(resType,theResType,NAMESIZE+1);
00857   strncpy(atomName,theAtomName,NAMESIZE+1);
00858   charge = theCharge;
00859   mass = theMass;
00860 }
00861 
00862 /* initializes this to be a bond between atoms <i> and <j> of type
00863    <type>, with <next> pointing to the next bond in the list. */
00864 GenericBond::GenericBond(int i, int j, int theType) {
00865   atomi=i;
00866   atomj=j;
00867   type=theType;
00868 }
00869 
00870 /* initializes this to be a angle between atoms <i>, <j>, and <k> of
00871    type <type>, with <next> pointing to the next angle in the list. */
00872 GenericAngle::GenericAngle(int i, int j, int k, int theType) {
00873   atomi=i;
00874   atomj=j;
00875   atomk=k;
00876   type=theType;
00877 }
00878 
00879 /* initializes this to be a angle between atoms <i>, <j>, <k>, and
00880    <l> of type <type> */
00881 GenericDihedral::GenericDihedral(int i, int j, int k, int l, int theType) {
00882   atomi=i;
00883   atomj=j;
00884   atomk=k;
00885   atoml=l;
00886   type=theType;
00887 }
00888 
00889 /* adds a bond to the list */
00890 void GenericMol::addBond(int atomi, int atomj, int type) {
00891   bondList.add(new GenericBond(atomi, atomj, type));
00892 }
00893 
00894 /* adds an angle to the list */
00895 void GenericMol::addAngle(int atomi, int atomj, int atomk, int type) {
00896   angleList.add(new GenericAngle(atomi, atomj, atomk, type));
00897 }
00898 
00899 /* adds a dihedral to the list */
00900 void GenericMol::addDihedral(int atomi, int atomj, int atomk,
00901                              int atoml, int type) {
00902   dihedralList.add(new GenericDihedral(atomi, atomj, atomk,
00903                                              atoml, type));
00904 }
00905 
00906 /* adds an atom to the list */
00907 void GenericMol::addAtom(const char *theType, int theTypeNum, int theResNum,
00908                          const char *theResType,const char *theAtomName,
00909                          Real theCharge, Real theMass) {
00910 
00911   atomList.add(new GenericAtom(theType,theTypeNum,theResNum,theResType,
00912                            theAtomName,theCharge,theMass));
00913 }
00914 
00915 /* initializes this to be the molecule with name <name> */
00916 GenericMol::GenericMol(const char *theName) {
00917   name = strdup(theName);
00918 }
00919 
00920 /* gets a bond */
00921 const GenericBond *GenericMol::getBond(int n) const {
00922   /* double-check */
00923   if(n >= bondList.size() || n<0) {
00924     fprintf(stderr,"Bond index %d out of bounds.\n",n);
00925     exit(1);
00926   }
00927   return bondList[n];
00928 }
00929 
00930 /* gets a angle */
00931 const GenericAngle *GenericMol::getAngle(int n) const {
00932   /* double-check */
00933   if(n >= angleList.size() || n<0) {
00934     fprintf(stderr,"Angle index %d out of bounds.\n",n);
00935     exit(1);
00936   }
00937   return angleList[n];
00938 }
00939 
00940 /* gets a dihedral */
00941 const GenericDihedral *GenericMol::getDihedral(int n) const {
00942   /* double-check */
00943   if(n >= dihedralList.size() || n<0) {
00944     fprintf(stderr,"Dihedral index %d out of bounds.\n",n);
00945     exit(1);
00946   }
00947   return dihedralList[n];
00948 }
00949 
00950 /* gets an atom */
00951 const GenericAtom *GenericMol::getAtom(int n) const {
00952   /* double-check */
00953   if(n >= atomList.size() || n<0) {
00954     fprintf(stderr,"Atom index %d out of bounds for %s.\n",n,name);
00955     exit(1);
00956   }
00957   return atomList[n];
00958 }
00959 
00960 /* initializes this to represent <theNum> copies of <theMol> */
00961 MolInst::MolInst(const GenericMol *theMol, int theNum) {
00962   mol = theMol;
00963   num = theNum;
00964 }
00965 
00966 /* get the total number of various things here */
00967 int MolInst::getNumAtoms() const {
00968   return getMol()->getNumAtoms() * getNum();
00969 }
00970 int MolInst::getNumBonds() const {
00971   return getMol()->getNumBonds() * getNum();
00972 }
00973 int MolInst::getNumAngles() const {
00974   return getMol()->getNumAngles() * getNum();
00975 }
00976 int MolInst::getNumDihedrals() const {
00977   return getMol()->getNumDihedrals() * getNum();
00978 }
00979 int MolInst::getNumRes() const {
00980   return getMol()->getNumRes()
00981     * getNum();
00982 }
00983 
00984 /* this gets the index of a bond in the table (adding an entry if
00985    none exists).
00986    b0: the natural bond length in A.
00987    kB: the spring constant in kcal/mol/A^2, where E=kx^2 to 1st order.
00988    funct: 1 for harmonic potential, 2 for fourth-order GROMACS96
00989    potential. */
00990 int BondTable::getIndex(float b0, float kB, int funct) {
00991   /* check to see if it is in the table already */
00992   int i;
00993   for(i=0;i<b0Array.size();i++) {
00994     if(fabs(b0-b0Array[i])<0.00001 &&
00995        fabs(kB-kBArray[i])<0.00001 &&
00996        funct == functArray[i]) {
00997       return i;
00998     }
00999   }
01000   
01001   /* nope, it wasn't in the table add a new element! */
01002   b0Array.add(b0);
01003   kBArray.add(kB);
01004   functArray.add(funct);
01005   typeaArray.add(NULL);
01006   typebArray.add(NULL);
01007   return b0Array.size()-1;
01008 }
01009 
01010 /* this gets the index of a angle in the table (adding an entry if
01011    none exists).
01012    b0: the natural angle length in A.
01013    kB: the spring constant in kcal/mol/A^2, where E=kx^2 to 1st order.
01014    funct: 1 for harmonic potential, 2 for fourth-order GROMACS96
01015    potential. */
01016 int AngleTable::getIndex(float th0, float kth, int funct) {
01017   /* check to see if it is in the table already */
01018   int i;
01019   for(i=0;i<th0Array.size();i++) {
01020     if(fabs(th0-th0Array[i])<0.00001 &&
01021        fabs(kth-kthArray[i])<0.00001 &&
01022        funct == functArray[i]) {
01023       return i;
01024     }
01025   }
01026   
01027   /* nope, it wasn't in the table add a new element! */
01028   th0Array.add(th0);
01029   kthArray.add(kth);
01030   functArray.add(funct);
01031   typeaArray.add(NULL);
01032   typebArray.add(NULL);
01033   typecArray.add(NULL);
01034   return th0Array.size()-1;
01035 }
01036 
01037 /* this function gets the index of a dihedral angle in the table
01038    (adding an entry if none exists).  The required number of
01039    parameters (see notes on class DihedralTable) must be stored in
01040    the array <c> */
01041 int DihedralTable::getIndex(const float *c, int mult, int funct) {
01042   /* check to see if it is in the table already */
01043   int i,j,jmax;
01044   if(funct==1 || funct==2) { /* for these we only need two params */
01045     jmax=2;
01046   } else { /* for RB we need all 6 */
01047     jmax=6;
01048   }
01049 
01050   for(i=0;i<cArray.size();i++) {
01051     int mismatch=0;
01052     if(mult != multArray[i]) continue;
01053     if(funct != functArray[i]) continue;
01054 
01055     for(j=0;j<jmax;j++) {
01056       if(fabs(c[j]-cArray[i][j])>=0.00001) {
01057         mismatch=1;
01058         break;
01059       }
01060     }
01061     if(!mismatch) {
01062       /* all of the parameters matched */
01063       return i;
01064     }
01065   }
01066   
01067   /* nope, it wasn't in the table add a new element! */
01068   addType(NULL,NULL,c,mult,funct);
01069   return cArray.size()-1;
01070 }
01071 
01072 /* getBondParams puts the parameters for bond-type <num> into the
01073    spaces pointed to by the other arguments. 
01074    b0 - natural length in A
01075    kB - spring constant in kcal/A^2 - E=kx^2
01076    funct - 1 for harmonic, 2 for special fourth-order GROMOS96 thing */
01077 void GromacsTopFile::getBondParams(int num, float *b0, float *kB, int *funct)
01078   const {
01079   bondTable.getParams(num,b0,kB,funct);
01080 }
01081 
01082 /* getAngleParams puts the parameters for angle-type <num> into the
01083    spaces pointed to by the other arguments.  */
01084 void GromacsTopFile::getAngleParams(int num, float *th0, float *kth, int *funct)
01085   const {
01086   angleTable.getParams(num,th0,kth,funct);
01087 }
01088 
01089 void GromacsTopFile::getDihedralParams(int num, float *c, int *mult, int *funct)
01090   const {
01091   dihedralTable.getParams(num,c,mult,funct);
01092 }
01093 
01094 void GromacsTopFile::getPairLJArrays2(int *indexA, int *indexB, Real *pairC6, Real *pairC12) {
01095   pairTable.getPairLJArrays2(indexA, indexB, pairC6, pairC12);
01096 }
01097 
01098 void GromacsTopFile::getPairGaussArrays2(int *indexA, int *indexB, Real *gaussA, Real *gaussMu1,
01099                                     Real *gaussSigma1, Real *gaussMu2, Real *gaussSigma2,
01100                                     Real *gaussRepulsive) {
01101   pairTable.getPairGaussArrays2(indexA, indexB, gaussA, gaussMu1, gaussSigma1, gaussMu2, gaussSigma2, gaussRepulsive);
01102 }
01103 
01104 /* getParams puts the parameters for bond-type <num> into the
01105    spaces pointed to by the other arguments. 
01106    b0 - natural length in A
01107    kB - spring constant in kcal/A^2 - E=kx^2
01108    funct - 1 for harmonic, 2 for special fourth-order GROMOS96 thing */
01109 void BondTable::getParams(int num, float *b0, float *kB, int *funct)
01110   const {
01111   *b0=b0Array[num];
01112   *kB=kBArray[num];
01113   *funct=functArray[num];
01114 }
01115 
01116 /* getParams puts the parameters for angle-type <num> into the
01117    spaces pointed to by the other arguments.  */
01118 void AngleTable::getParams(int num, float *th0, float *kth, int *funct)
01119   const {
01120   *th0=th0Array[num];
01121   *kth=kthArray[num];
01122   *funct=functArray[num];
01123 }
01124 
01125 /* getParams puts the parameters for angle-type <num> into the
01126    spaces pointed to by the other arguments.  The required number of
01127    parameters (see notes on class DihedralTable) will be stored in
01128    the array <c>, so make sure that c has size >= 6 */
01129 void DihedralTable::getParams(int num, float *c, int *mult, int *funct) const {
01130   int i;
01131 
01132   *funct=functArray[num]; /* first set the function */
01133 
01134   if(*funct==1 || *funct==2) { /* for these we only need two params */
01135     c[0]=cArray[num][0];
01136     c[1]=cArray[num][1];
01137   } else if(*funct==3) { /* for RB we need all 6 */
01138     for(i=0;i<6;i++) {
01139       c[i]=cArray[num][i];
01140     }
01141   } else { /* error */
01142     fprintf(stderr,"Bad function number %d - don't know what to do!\n",
01143             *funct);
01144     exit(1);
01145   }
01146 
01147   if(*funct==1) { /* return the multiplicity */
01148     *mult=multArray[num];
01149   }
01150 }
01151 
01152 /* this adds a entry for a angle type to the table, including the
01153    three atoms involved in the angle */
01154 void AngleTable::addType(const char *typea, const char *typeb,
01155                   const char *typec, float th0, float kth, int funct) {
01156   typeaArray.add(strdup(typea));
01157   typebArray.add(strdup(typeb));
01158   typecArray.add(strdup(typec));
01159   th0Array.add(th0);
01160   kthArray.add(kth);
01161   functArray.add(funct);
01162 }
01163 
01164 /* this adds a entry for a angle type to the table, including the
01165    two atoms involved in the angle.  The required number of
01166    parameters (see notes on class DihedralTable) must be stored in
01167    the array <c>.  Note that these two angles really refer to either
01168    atoms A and D or B and C depending on the dihedral type.  */
01169 void DihedralTable::addType(const char *typea, const char *typeb,
01170                          const float *c, int mult, int funct) {
01171   float *cNew;
01172   int i;
01173 
01174   if(typea != NULL) typeaArray.add(strdup(typea));
01175   if(typeb != NULL) typebArray.add(strdup(typeb));
01176   functArray.add(funct);
01177 
01178   if(funct==1) multArray.add(mult);
01179   else multArray.add(0);
01180 
01181   if(funct==1 || funct==2) { /* for these we only need two params */
01182     cNew = new float[2];
01183     cNew[0]=c[0];
01184     cNew[1]=c[1];
01185   } else if(funct==3) { /* for RB we need all 6 */
01186     cNew = new float[6];
01187     for(i=0;i<6;i++) {
01188       cNew[i]=c[i];
01189     }
01190   } else { /* error */
01191     fprintf(stderr,"Bad function number %d - don't know what to do!\n",
01192             funct);
01193     exit(1);
01194   }
01195   
01196   /* now push the actual parameters */
01197   cArray.add(cNew);
01198 }
01199 
01200 /* This version looks up the angle by atom type - the direction of
01201    the types doesn't matter (it can be A--B--C or C--B--A)!  If the
01202    specified angle/function combination cannot be found, it returns
01203    -1, otherwise it returns the index of the angletype. */
01204 int AngleTable::getParams(const char *typea, const char *typeb,
01205                    const char *typec, int funct,
01206                    float *th0, float *kth) const {
01207   int i;
01208   for(i=0;i<th0Array.size();i++) {
01209     if(typeaArray[i] == NULL || typebArray[i] == NULL
01210        || typecArray[i] == NULL) continue;
01211     if( (0==strcmp(typea,typeaArray[i]) &&  /* A--B--C */
01212          0==strcmp(typeb,typebArray[i]) &&
01213          0==strcmp(typec,typecArray[i]))    /* or */
01214      || (0==strcmp(typec,typeaArray[i]) &&
01215          0==strcmp(typeb,typebArray[i]) &&  /* C--B--A */
01216          0==strcmp(typea,typecArray[i]))
01217      && funct==functArray[i]) {
01218       *th0 = th0Array[i];
01219       *kth = kthArray[i];
01220       return i;
01221     }
01222   }
01223   return -1;
01224 }
01225 
01226 /* see the (long) notes in the prototype definition */
01227 int DihedralTable::getParams(const char *typea, const char *typeb,
01228                              const char *typec, const char *typed,
01229                              int funct, float *c, int *mult) const {
01230   int i,j;
01231   const char *comparea, *compareb;
01232   
01233   if(funct == 1 || funct == 3) { /* for these, we use the inner atoms */
01234     comparea = typeb;
01235     compareb = typec;
01236   } else { /* use the outer atoms */
01237     comparea = typea;
01238     compareb = typed;
01239   }
01240 
01241   for(i=0;i<cArray.size();i++) {
01242     if(typeaArray[i] == NULL || typebArray[i] == NULL)
01243       continue; /* no atom types specified */
01244     if(functArray[i] != funct)
01245       continue; /* wrong function type */
01246 
01247     if( (0==strcmp(comparea,typeaArray[i]) &&  /* A--B */
01248          0==strcmp(compareb,typebArray[i]))    /*  or  */
01249      || (0==strcmp(compareb,typeaArray[i]) &&
01250          0==strcmp(comparea,typebArray[i]))    /* B--A */
01251           ) {
01252       if(funct==1 || funct==2) { /* for these we only need two params */
01253         c[0]=cArray[i][0];
01254         c[1]=cArray[i][1];
01255         if(funct==1) {
01256           *mult = multArray[i];
01257         }
01258       } else if(funct==3) { /* for RB we need all 6 */
01259         for(j=0;j<6;j++) {
01260           c[j]=cArray[i][j];
01261         }
01262       }
01263       return i;
01264     }
01265   }
01266   return -1;
01267 }
01268 
01269 /* this adds a entry for a bond type to the table, including the two
01270    atoms involved in the bond */
01271 void BondTable::addType(const char *typea, const char *typeb,
01272                             float b0, float kB, int funct) {
01273   b0Array.add(b0);
01274   kBArray.add(kB);
01275   functArray.add(funct);
01276   typeaArray.add(strdup(typea));
01277   typebArray.add(strdup(typeb));
01278 }
01279 
01280 /* This version looks up the bond by atom type - the order of the
01281    types doesn't matter! */
01282 int BondTable::getParams(const char *typea, const char *typeb,
01283                               int funct, float *b0, float *kB) const {
01284   int i;
01285   for(i=0;i<b0Array.size();i++) {
01286     if(typeaArray[i] == NULL || typebArray[i] == NULL) continue;
01287     if( (0==strcmp(typea,typeaArray[i]) &&
01288          0==strcmp(typeb,typebArray[i]))
01289      || (0==strcmp(typeb,typeaArray[i]) &&
01290          0==strcmp(typea,typebArray[i]))
01291      && funct==functArray[i]) {
01292       *b0 = b0Array[i];
01293       *kB = kBArray[i];
01294       return i;
01295     }
01296   }
01297   return -1;
01298 }
01299 
01300 /* this adds a entry for an atom type to the table */
01301 void AtomTable::addType(const char *type, float m, float q,
01302                             float c6, float c12) {
01303   typeArray.add(strdup(type));
01304   mArray.add(m);
01305   qArray.add(q);
01306   c6Array.add(c6);
01307   c12Array.add(c12);
01308 }
01309 
01310 /* This looks up the atom type by number, returning it in the string
01311    <type>, which must be at least 11 characters long. */
01312 void AtomTable::getType(int num, char *type) const {
01313   if(num>=mArray.size() || num<0) {
01314     fprintf(stderr,"atomParam index %d out of bounds!\n",num);
01315     exit(1);
01316   }
01317   strncpy(type,typeArray[num],NAMESIZE+1);
01318 }
01319 
01320 /* This looks up the atom by type - if it cannot be found, this
01321    returns -1, otherwise this returns the index of the atomtype (for
01322    consistency - this is really a useless number.) */
01323 int AtomTable::getParams(const char *type, float *m, float *q,
01324                   float *c6, float *c12) const {
01325   int i;
01326   for(i=0;i<mArray.size();i++) {
01327     if(typeArray[i]==NULL) {
01328       fprintf(stderr,"Found NULL atom type in list.\n");
01329       exit(1);
01330     }
01331     if(0==strcmp(typeArray[i],type)) {
01332       *m = mArray[i];
01333       *q = qArray[i];
01334       *c6 = c6Array[i];
01335       *c12 = c12Array[i];
01336       return i;
01337     }
01338   }
01339   return -1;
01340 }
01341 
01342 /* finds the index from the two interacting types, or returns -1 */
01343 int VDWTable::getIndex(const char *typea, const char *typeb) const {
01344   int i;
01345   for(i=0;i<c6Array.size();i++) {
01346     if((0==strcmp(typea,typeAArray[i]) &&
01347         0==strcmp(typeb,typeBArray[i])) ||
01348        (0==strcmp(typeb,typeAArray[i]) &&
01349         0==strcmp(typea,typeBArray[i]))) {
01350       return i;
01351     }
01352   }
01353   return -1;
01354 }
01355 
01356 /* these add VDW parameters to the list */
01357 void VDWTable::addType(const char *typea, const char *typeb, Real c6, 
01358                        Real c12) {
01359   int i;
01360 
01361   /* check to see if the pair is already in the table */
01362   i = getIndex(typea,typeb);
01363   if(i != -1) {  /* it was in the table */
01364     c6Array[i] = c6;
01365     c12Array[i] = c12;
01366   }
01367   else { /* it wasn't in the table - add it! */
01368     typeAArray.add(strdup(typea));
01369     typeBArray.add(strdup(typeb));
01370     c6Array.add(c6);
01371     c12Array.add(c12);
01372     c6PairArray.add(0);
01373     c12PairArray.add(0);
01374   }
01375 }
01376 
01377 void VDWTable::add14Type(const char *typea, const char *typeb,
01378                Real c6pair, Real c12pair) {
01379   int i;
01380 
01381   /* check to see if the pair is already in the table */
01382   i = getIndex(typea,typeb);
01383   if(i != -1) {  /* it was in the table */
01384     c6PairArray[i] = c6pair;
01385     c12PairArray[i] = c12pair;
01386   }
01387   else { /* it wasn't in the table - add it! */
01388     typeAArray.add(strdup(typea));
01389     typeBArray.add(strdup(typeb));
01390     c6PairArray.add(c6pair);
01391     c12PairArray.add(c12pair);
01392     c6Array.add(0);
01393     c12Array.add(0);
01394   }
01395 }
01396 
01397 /* this returns the VDW interaction parameters that were added to
01398    the list earlier, and returns the index or the parameters (a
01399    useless number) or -1 if it can't find the specified types. */
01400 int VDWTable::getParams(const char *typea, const char *typeb,
01401               Real *c6, Real *c12, Real *c6pair, Real *c12pair) const {
01402   int i;
01403   /* check to see if the pair is already in the table */
01404   i = getIndex(typea,typeb);
01405   if(i != -1) {  /* it was in the table - return the parameters  */
01406     *c6 = c6Array[i];
01407     *c12 = c12Array[i];
01408     *c6pair = c6PairArray[i];
01409     *c12pair = c12PairArray[i];
01410   }
01411   return i;
01412 }
01413 
01414 /* getVDWParams returs the Lennard-Jones bonding parameters for the
01415    specified two atom types, and the modified bonding parameters
01416    for 1-4 L-J interactons (c6pair, c12pair). */
01417 void GromacsTopFile::getVDWParams(int numa, int numb,
01418                   Real *c6, Real *c12, Real *c6pair, Real *c12pair) const {
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 }
01467 
01468 /* Start of JLai modifications August 16th, 2012 */
01469 int PairTable::addPairLJType2(int indexA, int indexB, Real pairC6, Real pairC12) {
01470   GroLJPair glp;
01471   glp.indxLJA = indexA;
01472   glp.indxLJB = indexB;
01473   glp.c6pair = pairC6;
01474   glp.c12pair = pairC12;
01475   pairlistLJ.add(glp);
01476   numLJPair++;
01477   return 0;
01478 
01479   // Insert the second copy
01480   GroLJPair glp2;
01481   glp2.indxLJA = indexB;
01482   glp2.indxLJB = indexA;
01483   glp2.c6pair = pairC6;
01484   glp2.c12pair = pairC12;
01485   pairlistLJ.add(glp2);
01486   numLJPair++;
01487   return 0;
01488 }
01489 
01490 int PairTable::addPairGaussType2(int indexA, int indexB, Real gaussA, Real gaussMu1,
01491                                 Real gaussSigma1) {
01492   return addPairGaussType2(indexA, indexB, gaussA, gaussMu1, gaussSigma1, 0.0, 0.0, 0.0);
01493 }
01494 
01495 int PairTable::addPairGaussType2(int indexA, int indexB, Real gaussA, Real gaussMu1,
01496                              Real gaussSigma1, Real gaussRepulsive) {
01497   return addPairGaussType2(indexA, indexB, gaussA, gaussMu1, gaussSigma1, 0.0, 0.0, gaussRepulsive);
01498 }
01499 
01500 int PairTable::addPairGaussType2(int indexA, int indexB, Real gaussA, Real gaussMu1,
01501                              Real gaussSigma1, Real gaussMu2, Real gaussSigma2, 
01502                              Real gaussRepulsive) {
01503   GroGaussPair ggp;
01504   ggp.indxGaussA = indexA;
01505   ggp.indxGaussB = indexB;
01506   ggp.gA = gaussA;
01507   ggp.gMu1 = gaussMu1;
01508   ggp.giSigma1 = gaussSigma1;
01509   ggp.gMu2 = 0.0;
01510   ggp.giSigma2 = 0.0;
01511   ggp.gRepulsive = 0.0;
01512   pairlistGauss.add(ggp);
01513   numGaussPair++;
01514   GroGaussPair ggp2;
01515   ggp2.indxGaussA = indexB;
01516   ggp2.indxGaussB = indexA;
01517   ggp2.gA = gaussA;
01518   ggp2.gMu1 = gaussMu1;
01519   ggp2.giSigma1 = gaussSigma1;
01520   ggp2.gMu2 = 0.0;
01521   ggp2.giSigma2 = 0.0;
01522   ggp2.gRepulsive = 0.0;
01523   numGaussPair++;
01524   pairlistGauss.add(ggp2);
01525   return 0;
01526 }
01527 
01528 void PairTable::getPairLJArrays2(int *indexA, int *indexB, Real *pairC6, Real *pairC12) {
01529 
01530   std::sort(pairlistLJ.begin(),pairlistLJ.end(),GroLJCompare);
01531   ResizeArray<GroLJPair>::iterator it;
01532   for(int i = 0; i < numLJPair; i++) {  
01533     indexA[i] = pairlistLJ[i].indxLJA;
01534     indexB[i] = pairlistLJ[i].indxLJB;
01535     pairC6[i] = pairlistLJ[i].c6pair;
01536     pairC12[i] = pairlistLJ[i].c12pair;
01537     }
01538 }
01539 
01540 void PairTable::getPairGaussArrays2(int *indexA, int *indexB, Real *gaussA, Real *gaussMu1,
01541                                    Real *gaussSigma1, Real *gaussMu2, Real *gaussSigma2,
01542                                    Real *gaussRepulsive) {
01543   std::sort(pairlistGauss.begin(),pairlistGauss.end(),GroGaussCompare);
01544   for(int i = 0; i < numGaussPair; i++) {
01545     indexA[i] = pairlistGauss[i].indxGaussA;
01546     indexB[i] = pairlistGauss[i].indxGaussB;
01547     gaussA[i] = pairlistGauss[i].gA;
01548     gaussMu1[i] = pairlistGauss[i].gMu1;
01549     gaussSigma1[i] = pairlistGauss[i].giSigma1;
01550     gaussMu2[i] = pairlistGauss[i].gMu2;
01551     gaussSigma2[i] = pairlistGauss[i].giSigma2;
01552     gaussRepulsive[i] = pairlistGauss[i].gRepulsive;
01553   }
01554 }
01555 
01556 bool PairTable::GroLJCompare (GroLJPair A, GroLJPair B) {
01557   if(A.indxLJA < B.indxLJA) {
01558     return true;
01559   } else if(A.indxLJA == B.indxLJA) {
01560     return (A.indxLJB < B.indxLJB);
01561   } 
01562   return false;
01563 }
01564 
01565 bool PairTable::GroGaussCompare (GroGaussPair A, GroGaussPair B) {
01566   if(A.indxGaussA < B.indxGaussA) {
01567     return true;
01568   } else if(A.indxGaussA == B.indxGaussA) {
01569     return (A.indxGaussB < B.indxGaussB);
01570   } 
01571   return false;
01572 }
01573 /* End of JLai Modifications */

Generated on Sat Sep 23 01:17:13 2017 for NAMD by  doxygen 1.4.7