Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

CompressPsf.C

Go to the documentation of this file.
00001 #include <algorithm>
00002 #include "CompressPsf.h"
00003 #include "strlib.h"
00004 #include "Molecule.h"
00005 #include "Parameters.h"
00006 #include "SimParameters.h"
00007 #include "InfoStream.h"
00008 #include "UniqueSet.h"
00009 #include "UniqueSetIter.h"
00010 
00011 
00030 //global variables recording compressed psf file information
00031 Molecule *g_mol = NULL;
00032 Parameters *g_param = NULL;
00033 SimParameters *g_simParam = NULL;
00034 ConfigList *g_cfgList = NULL;
00035 
00036 struct BasicAtomInfo
00037 {
00038     short segNameIdx;
00039     short resNameIdx;
00040     short atomNameIdx;
00041     short atomTypeIdx;
00042     short chargeIdx;
00043     short massIdx;
00044     int resID;
00045 
00046     int atomSigIdx;
00047     int exclSigIdx;
00048 };
00049 
00050 struct AtomSigInfo
00051 {
00052     vector<short> bondSigIndices;
00053     vector<short> angleSigIndices;
00054     vector<short> dihedralSigIndices;
00055     vector<short> improperSigIndices;
00056     vector<short> crosstermSigIndices;
00057 
00058     AtomSigInfo()
00059     {}
00060     AtomSigInfo(const AtomSigInfo& sig)
00061     {
00062         bondSigIndices.clear();
00063         for(int i=0; i<sig.bondSigIndices.size(); i++)
00064             bondSigIndices.push_back(sig.bondSigIndices[i]);
00065 
00066         angleSigIndices.clear();
00067         for(int i=0; i<sig.angleSigIndices.size(); i++)
00068             angleSigIndices.push_back(sig.angleSigIndices[i]);
00069 
00070         dihedralSigIndices.clear();
00071         for(int i=0; i<sig.dihedralSigIndices.size(); i++)
00072             dihedralSigIndices.push_back(sig.dihedralSigIndices[i]);
00073 
00074         improperSigIndices.clear();
00075         for(int i=0; i<sig.improperSigIndices.size(); i++)
00076             improperSigIndices.push_back(sig.improperSigIndices[i]);
00077 
00078         crosstermSigIndices.clear();
00079         for(int i=0; i<sig.crosstermSigIndices.size(); i++)
00080             crosstermSigIndices.push_back(sig.crosstermSigIndices[i]);
00081     }
00082 
00083     ~AtomSigInfo()
00084     {
00085         bondSigIndices.clear();
00086         angleSigIndices.clear();
00087         dihedralSigIndices.clear();
00088         improperSigIndices.clear();
00089         crosstermSigIndices.clear();
00090     }
00091 
00092     void sortTupleSigIndices()
00093     {
00094         sort(bondSigIndices.begin(), bondSigIndices.end());
00095         sort(angleSigIndices.begin(), angleSigIndices.end());
00096         sort(dihedralSigIndices.begin(), dihedralSigIndices.end());
00097         sort(improperSigIndices.begin(), improperSigIndices.end());
00098         sort(crosstermSigIndices.begin(), crosstermSigIndices.end());
00099     }
00100 };
00101 
00102 int operator==(const AtomSigInfo &s1, const AtomSigInfo& s2)
00103 {
00104     if(s1.bondSigIndices.size() != s2.bondSigIndices.size())
00105         return 0;
00106     if(s1.angleSigIndices.size() != s2.angleSigIndices.size())
00107         return 0;
00108     if(s1.dihedralSigIndices.size() != s2.dihedralSigIndices.size())
00109         return 0;
00110     if(s1.improperSigIndices.size() != s2.improperSigIndices.size())
00111         return 0;
00112     if(s1.crosstermSigIndices.size() != s2.crosstermSigIndices.size())
00113         return 0;
00114 
00115     int equalCnt;
00116     equalCnt=0;
00117     int bondSigCnt = s1.bondSigIndices.size();
00118     for(int i=0; i<bondSigCnt; i++)
00119         equalCnt += (s1.bondSigIndices[i]==s2.bondSigIndices[i]);
00120     if(equalCnt!=bondSigCnt)
00121         return 0;
00122 
00123     equalCnt=0;
00124     int angleSigCnt = s1.angleSigIndices.size();
00125     for(int i=0; i<angleSigCnt; i++)
00126         equalCnt += (s1.angleSigIndices[i]==s2.angleSigIndices[i]);
00127     if(equalCnt!=angleSigCnt)
00128         return 0;
00129 
00130     equalCnt=0;
00131     int dihedralSigCnt = s1.dihedralSigIndices.size();
00132     for(int i=0; i<dihedralSigCnt; i++)
00133         equalCnt += (s1.dihedralSigIndices[i]==s2.dihedralSigIndices[i]);
00134     if(equalCnt!=dihedralSigCnt)
00135         return 0;
00136 
00137     equalCnt=0;
00138     int improperSigCnt = s1.improperSigIndices.size();
00139     for(int i=0; i<improperSigCnt; i++)
00140         equalCnt += (s1.improperSigIndices[i]==s2.improperSigIndices[i]);
00141     if(equalCnt!=improperSigCnt)
00142         return 0;
00143 
00144     equalCnt=0;
00145     int crosstermSigCnt = s1.crosstermSigIndices.size();
00146     for(int i=0; i<crosstermSigCnt; i++)
00147         equalCnt += (s1.crosstermSigIndices[i]==s2.crosstermSigIndices[i]);
00148     if(equalCnt!=crosstermSigCnt)
00149         return 0;
00150 
00151     return 1;
00152 }
00153 
00154 struct ExclSigInfo
00155 {
00156     vector<int> fullExclOffset;
00157     vector<int> modExclOffset;
00158 
00159     ExclSigInfo()
00160     {}
00161     ExclSigInfo(const ExclSigInfo& sig)
00162     {
00163         fullExclOffset.clear();
00164         for(int i=0; i<sig.fullExclOffset.size(); i++)
00165             fullExclOffset.push_back(sig.fullExclOffset[i]);
00166 
00167         modExclOffset.clear();
00168         for(int i=0; i<sig.modExclOffset.size(); i++)
00169             modExclOffset.push_back(sig.modExclOffset[i]);
00170     }
00171 
00172     ~ExclSigInfo()
00173     {
00174         fullExclOffset.clear();
00175         modExclOffset.clear();
00176     }
00177 
00178     void sortExclOffset()
00179     {
00180         sort(fullExclOffset.begin(), fullExclOffset.end());
00181         sort(modExclOffset.begin(), modExclOffset.end());
00182     }
00183 };
00184 int operator==(const ExclSigInfo &s1, const ExclSigInfo &s2)
00185 {
00186     if(s1.fullExclOffset.size()!=s2.fullExclOffset.size())
00187         return 0;
00188     if(s1.modExclOffset.size()!=s2.modExclOffset.size())
00189         return 0;
00190 
00191     for(int i=0; i<s1.fullExclOffset.size(); i++)
00192     {
00193         if(s1.fullExclOffset[i] != s2.fullExclOffset[i])
00194             return 0;
00195     }
00196 
00197     for(int i=0; i<s1.modExclOffset.size(); i++)
00198     {
00199         if(s1.modExclOffset[i] != s2.modExclOffset[i])
00200             return 0;
00201     }
00202     return 1;
00203 }
00204 
00205 
00206 vector<string> segNamePool;
00207 vector<string> resNamePool;
00208 vector<string> atomNamePool;
00209 vector<string> atomTypePool;
00210 vector<Real> chargePool;
00211 vector<Real> massPool;
00212 vector<AtomSigInfo> atomSigPool;
00213 BasicAtomInfo *atomData;
00214 
00215 //Recording cluster information after reading all bonds info
00216 int *eachAtomClusterID = NULL;
00217 vector<int> eachClusterSize;
00218 
00219 vector<TupleSignature> sigsOfBonds;
00220 vector<TupleSignature> sigsOfAngles;
00221 vector<TupleSignature> sigsOfDihedrals;
00222 vector<TupleSignature> sigsOfImpropers;
00223 vector<TupleSignature> sigsOfCrossterms;
00224 AtomSigInfo *eachAtomSigs;
00225 
00226 vector<ExclSigInfo> sigsOfExclusions;
00227 ExclSigInfo *eachAtomExclSigs;
00228 
00229 //Structures for extraBond options
00230 vector<Bond> extraBonds;
00231 vector<Angle> extraAngles;
00232 vector<Dihedral> extraDihedrals;
00233 vector<Improper> extraImpropers;
00234 vector<BondValue> extraBondParams;
00235 vector<AngleValue> extraAngleParams;
00236 vector<DihedralValue> extraDihedralParams;
00237 vector<ImproperValue> extraImproperParams;
00238 
00239 int operator==(const BondValue &b1, const BondValue &b2)
00240 {
00241     return (b1.k==b2.k) && (b1.x0==b2.x0);
00242 }
00243 
00244 int operator==(const AngleValue &a1, const AngleValue &a2)
00245 {
00246     return (a1.k==a2.k) && (a1.k_ub==a2.k_ub) && (a1.r_ub==a2.r_ub) && (a1.theta0==a2.theta0);
00247 }
00248 
00249 int operator!=(const FourBodyConsts& f1, const FourBodyConsts& f2)
00250 {
00251     return (f1.delta!=f2.delta) || (f1.k!=f2.k) || (f1.n!=f2.n);
00252 }
00253 
00254 int operator==(const DihedralValue &d1, const DihedralValue &d2)
00255 {
00256     if(d1.multiplicity != d2.multiplicity)
00257         return 0;
00258     for(int i=0; i<MAX_MULTIPLICITY; i++)
00259     {
00260         if(d1.values[i] != d2.values[i])
00261             return 0;
00262     }
00263     return 1;
00264 }
00265 
00266 int operator==(const ImproperValue &d1, const ImproperValue &d2)
00267 {
00268     if(d1.multiplicity != d2.multiplicity)
00269         return 0;
00270     for(int i=0; i<MAX_MULTIPLICITY; i++)
00271     {
00272         if(d1.values[i] != d2.values[i])
00273             return 0;
00274     }
00275     return 1;
00276 }
00277 
00278 void readPsfFile(char *psfFileName);
00279 void integrateAllAtomSigs();
00280 void outputPsfFile(FILE *ofp);
00281 
00282 //reading extraBond's information
00283 void getExtraBonds(StringList *file);
00284 
00285 //reading atom's basic information
00286 void getAtomData(FILE *ifp);
00287 void getBondData(FILE *ifp);
00288 void getAngleData(FILE *ifp);
00289 void getDihedralData(FILE *ifp);
00290 void getImproperData(FILE *ifp);
00291 void getDonorData(FILE *ifp);
00292 void getAcceptorData(FILE *ifp);
00293 void getExclusionData(FILE *ifp);
00294 void getCrosstermData(FILE *ifp);
00295 
00296 //Functions related with building exclusions
00297 void buildExclusions();
00298 void build12Excls(UniqueSet<Exclusion>&, vector<int> *);
00299 void build13Excls(UniqueSet<Exclusion>&, vector<int> *);
00300 void build14Excls(UniqueSet<Exclusion>&, vector<int> *, int);
00301 
00302 void clearGlobalVectors()
00303 {
00304     segNamePool.clear();
00305     resNamePool.clear();
00306     atomNamePool.clear();
00307     atomTypePool.clear();
00308     chargePool.clear();
00309     massPool.clear();
00310     sigsOfBonds.clear();
00311     sigsOfAngles.clear();
00312     sigsOfDihedrals.clear();
00313     sigsOfImpropers.clear();
00314     sigsOfExclusions.clear();
00315 
00316     eachClusterSize.clear();
00317 }
00318 
00319 void compress_psf_file(Molecule *mol, char *psfFileName, Parameters *param, SimParameters *simParam, ConfigList *cfgList)
00320 {
00321 
00322     g_mol = mol;
00323     g_param = param;
00324     g_simParam = simParam; //used for building exclusions
00325     g_cfgList = cfgList; //used for integrating extra bonds
00326 
00327     //read psf files
00328     readPsfFile(psfFileName);
00329 
00330     integrateAllAtomSigs();
00331 
00332     buildExclusions();
00333 
00334 
00335     char *outFileName = new char[strlen(psfFileName)+10];
00336     sprintf(outFileName, "%s.inter", psfFileName);
00337     FILE *ofp = fopen(outFileName, "w");
00338     delete [] outFileName;
00339     //output compressed psf file
00340     outputPsfFile(ofp);
00341     fclose(ofp);
00342 }
00343 
00344 //Almost same with Molecule::read_psf_file
00345 void readPsfFile(char *fname)
00346 {
00347     char err_msg[512];  //  Error message for NAMD_die
00348     char buffer[512];  //  Buffer for file reading
00349     int i;      //  Loop counter
00350     int NumTitle;    //  Number of Title lines in .psf file
00351     int ret_code;    //  ret_code from NAMD_read_line calls
00352     FILE *psf_file;
00353     Parameters *params = g_param;
00354 
00355     /* Try and open the .psf file           */
00356     if ( (psf_file = fopen(fname, "r")) == NULL)
00357     {
00358         sprintf(err_msg, "UNABLE TO OPEN .psf FILE %s", fname);
00359         NAMD_die(err_msg);
00360     }
00361 
00362     /*  Read till we have the first non-blank line of file    */
00363     ret_code = NAMD_read_line(psf_file, buffer);
00364 
00365     while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00366     {
00367         ret_code = NAMD_read_line(psf_file, buffer);
00368     }
00369 
00370     /*  Check to see if we dropped out of the loop because of a     */
00371     /*  read error.  This shouldn't happen unless the file is empty */
00372     if (ret_code!=0)
00373     {
00374         sprintf(err_msg, "EMPTY .psf FILE %s", fname);
00375         NAMD_die(err_msg);
00376     }
00377 
00378     /*  The first non-blank line should contain the word "psf".    */
00379     /*  If we can't find it, die.               */
00380     if (!NAMD_find_word(buffer, "psf"))
00381     {
00382         sprintf(err_msg, "UNABLE TO FIND \"PSF\" STRING IN PSF FILE %s",
00383                 fname);
00384         NAMD_die(err_msg);
00385     }
00386 
00387     /*  Read until we find the next non-blank line      */
00388     ret_code = NAMD_read_line(psf_file, buffer);
00389 
00390     while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00391     {
00392         ret_code = NAMD_read_line(psf_file, buffer);
00393     }
00394 
00395     /*  Check to see if we dropped out of the loop because of a     */
00396     /*  read error.  This shouldn't happen unless there is nothing  */
00397     /*  but the PSF line in the file        */
00398     if (ret_code!=0)
00399     {
00400         sprintf(err_msg, "MISSING EVERYTHING BUT PSF FROM %s", fname);
00401         NAMD_die(err_msg);
00402     }
00403 
00404     /*  This line should have the word "NTITLE" in it specifying    */
00405     /*  how many title lines there are        */
00406     if (!NAMD_find_word(buffer, "NTITLE"))
00407     {
00408         sprintf(err_msg,"CAN NOT FIND \"NTITLE\" STRING IN PSF FILE %s",
00409                 fname);
00410         NAMD_die(err_msg);
00411     }
00412 
00413     sscanf(buffer, "%d", &NumTitle);
00414 
00415     /*  Now skip the next NTITLE non-blank lines and then read in the*/
00416     /*  line which should contain NATOM        */
00417     i=0;
00418 
00419     while ( ((ret_code=NAMD_read_line(psf_file, buffer)) == 0) &&
00420             (i<NumTitle) )
00421     {
00422         if (!NAMD_blank_string(buffer))
00423             i++;
00424     }
00425 
00426     /*  Make sure we didn't exit because of a read error    */
00427     if (ret_code!=0)
00428     {
00429         sprintf(err_msg, "FOUND EOF INSTEAD OF NATOM IN PSF FILE %s",
00430                 fname);
00431         NAMD_die(err_msg);
00432     }
00433 
00434     while (NAMD_blank_string(buffer))
00435     {
00436         NAMD_read_line(psf_file, buffer);
00437     }
00438 
00439     /*  Check to make sure we have the line we want      */
00440     if (!NAMD_find_word(buffer, "NATOM"))
00441     {
00442         sprintf(err_msg, "DIDN'T FIND \"NATOM\" IN PSF FILE %s",
00443                 fname);
00444         NAMD_die(err_msg);
00445     }
00446 
00447     /*  Read in the number of atoms, and then the atoms themselves  */
00448     sscanf(buffer, "%d", &g_mol->numAtoms);
00449 
00450     getAtomData(psf_file);
00451 
00452     //read extra bonds/angles/dihedrals/impropers information first
00453     //and then integrate them into the following reading procedures
00454     if(g_simParam->extraBondsOn)
00455     {
00456         getExtraBonds(g_cfgList->find("extraBondsFile"));
00457     }
00458 
00459     //initialize eachAtomSigs
00460     eachAtomSigs = new AtomSigInfo[g_mol->numAtoms];
00461 
00462     /*  Read until we find the next non-blank line      */
00463     ret_code = NAMD_read_line(psf_file, buffer);
00464 
00465     while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00466     {
00467         ret_code = NAMD_read_line(psf_file, buffer);
00468     }
00469 
00470     /*  Check to make sure we didn't hit the EOF      */
00471     if (ret_code != 0)
00472     {
00473         NAMD_die("EOF ENCOUNTERED LOOKING FOR NBONDS IN PSF");
00474     }
00475 
00476     /*  Look for the string "NBOND"          */
00477     if (!NAMD_find_word(buffer, "NBOND"))
00478     {
00479         NAMD_die("DID NOT FIND NBOND AFTER ATOM LIST IN PSF");
00480     }
00481 
00482     /*  Read in the number of bonds and then the bonds themselves  */
00483     sscanf(buffer, "%d", &g_mol->numBonds);
00484 
00485     if (g_mol->numBonds)
00486         getBondData(psf_file);
00487 
00488     /*  Read until we find the next non-blank line      */
00489     ret_code = NAMD_read_line(psf_file, buffer);
00490 
00491     while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00492     {
00493         ret_code = NAMD_read_line(psf_file, buffer);
00494     }
00495 
00496     /*  Check to make sure we didn't hit the EOF      */
00497     if (ret_code != 0)
00498     {
00499         NAMD_die("EOF ENCOUNTERED LOOKING FOR NTHETA IN PSF");
00500     }
00501 
00502     /*  Look for the string "NTHETA"        */
00503     if (!NAMD_find_word(buffer, "NTHETA"))
00504     {
00505         NAMD_die("DID NOT FIND NTHETA AFTER BOND LIST IN PSF");
00506     }
00507 
00508     /*  Read in the number of angles and then the angles themselves */
00509     sscanf(buffer, "%d", &g_mol->numAngles);
00510 
00511     if (g_mol->numAngles)
00512         getAngleData(psf_file);
00513 
00514     /*  Read until we find the next non-blank line      */
00515     ret_code = NAMD_read_line(psf_file, buffer);
00516 
00517     while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00518     {
00519         ret_code = NAMD_read_line(psf_file, buffer);
00520     }
00521 
00522     /*  Check to make sure we didn't hit the EOF      */
00523     if (ret_code != 0)
00524     {
00525         NAMD_die("EOF ENCOUNTERED LOOKING FOR NPHI IN PSF");
00526     }
00527 
00528     /*  Look for the string "NPHI"          */
00529     if (!NAMD_find_word(buffer, "NPHI"))
00530     {
00531         NAMD_die("DID NOT FIND NPHI AFTER ANGLE LIST IN PSF");
00532     }
00533 
00534     /*  Read in the number of dihedrals and then the dihedrals      */
00535     sscanf(buffer, "%d", &g_mol->numDihedrals);
00536 
00537     if (g_mol->numDihedrals)
00538         getDihedralData(psf_file);
00539 
00540     /*  Read until we find the next non-blank line      */
00541     ret_code = NAMD_read_line(psf_file, buffer);
00542 
00543     while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00544     {
00545         ret_code = NAMD_read_line(psf_file, buffer);
00546     }
00547 
00548     /*  Check to make sure we didn't hit the EOF      */
00549     if (ret_code != 0)
00550     {
00551         NAMD_die("EOF ENCOUNTERED LOOKING FOR NIMPHI IN PSF");
00552     }
00553 
00554     /*  Look for the string "NIMPHI"        */
00555     if (!NAMD_find_word(buffer, "NIMPHI"))
00556     {
00557         NAMD_die("DID NOT FIND NIMPHI AFTER ATOM LIST IN PSF");
00558     }
00559 
00560     /*  Read in the number of Impropers and then the impropers  */
00561     sscanf(buffer, "%d", &g_mol->numImpropers);
00562 
00563     if (g_mol->numImpropers)
00564         getImproperData(psf_file);
00565 
00566     /*  Read until we find the next non-blank line      */
00567     ret_code = NAMD_read_line(psf_file, buffer);
00568 
00569     while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00570     {
00571         ret_code = NAMD_read_line(psf_file, buffer);
00572     }
00573 
00574     /*  Check to make sure we didn't hit the EOF      */
00575     if (ret_code != 0)
00576     {
00577         NAMD_die("EOF ENCOUNTERED LOOKING FOR NDON IN PSF");
00578     }
00579 
00580     /*  Look for the string "NDON"        */
00581     if (!NAMD_find_word(buffer, "NDON"))
00582     {
00583         NAMD_die("DID NOT FIND NDON AFTER ATOM LIST IN PSF");
00584     }
00585 
00586     /*  Read in the number of hydrogen bond donors and then the donors */
00587     sscanf(buffer, "%d", &g_mol->numDonors);
00588 
00589     if (g_mol->numDonors)
00590         getDonorData(psf_file);
00591 
00592     /*  Read until we find the next non-blank line      */
00593     ret_code = NAMD_read_line(psf_file, buffer);
00594 
00595     while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00596     {
00597         ret_code = NAMD_read_line(psf_file, buffer);
00598     }
00599 
00600     /*  Check to make sure we didn't hit the EOF      */
00601     if (ret_code != 0)
00602     {
00603         NAMD_die("EOF ENCOUNTERED LOOKING FOR NACC IN PSF");
00604     }
00605 
00606     /*  Look for the string "NACC"        */
00607     if (!NAMD_find_word(buffer, "NACC"))
00608     {
00609         NAMD_die("DID NOT FIND NACC AFTER ATOM LIST IN PSF");
00610     }
00611 
00612     /*  Read in the number of hydrogen bond donors and then the donors */
00613     sscanf(buffer, "%d", &g_mol->numAcceptors);
00614 
00615     if (g_mol->numAcceptors)
00616         getAcceptorData(psf_file);
00617 
00618     /*  look for the explicit non-bonded exclusion section.     */
00619     while (!NAMD_find_word(buffer, "NNB"))
00620     {
00621         ret_code = NAMD_read_line(psf_file, buffer);
00622 
00623         if (ret_code != 0)
00624         {
00625             NAMD_die("EOF ENCOUNTERED LOOKING FOR NNB IN PSF FILE");
00626         }
00627     }
00628 
00629     /*  Read in the number of exclusions and then the exclusions    */
00630     sscanf(buffer, "%d", &g_mol->numExclusions);
00631 
00632     if (g_mol->numExclusions)
00633         getExclusionData(psf_file);
00634 
00635     /*  look for the cross-term section.     */
00636     int crossterms_present = 1;
00637     while (!NAMD_find_word(buffer, "NCRTERM"))
00638     {
00639         ret_code = NAMD_read_line(psf_file, buffer);
00640 
00641         if (ret_code != 0)
00642         {
00643             // hit EOF before finding cross-term section
00644             crossterms_present = 0;
00645             break;
00646         }
00647     }
00648 
00649     if ( crossterms_present)
00650     {
00651 
00652         /*  Read in the number of cross-terms and then the cross-terms*/
00653         sscanf(buffer, "%d", &g_mol->numCrossterms);
00654 
00655         if (g_mol->numCrossterms)
00656             getCrosstermData(psf_file);
00657 
00658     }
00659 
00660     //  analyze the data and find the status of each atom
00661     //build_atom_status();
00662 
00663     fclose(psf_file);
00664 }
00665 
00666 void integrateAllAtomSigs()
00667 {
00668     printf("Bond sigs:  %d\n", (int)sigsOfBonds.size());
00669     printf("Angle sigs:  %d\n", (int)sigsOfAngles.size());
00670     printf("Dihedral sigs:  %d\n", (int)sigsOfDihedrals.size());
00671     printf("Improper sigs:  %d\n", (int)sigsOfImpropers.size());
00672     printf("Crossterm sigs:  %d\n", (int)sigsOfCrossterms.size());
00673 
00674 
00675     for(int i=0; i<g_mol->numAtoms; i++)
00676     {
00677         eachAtomSigs[i].sortTupleSigIndices();
00678         int poolIndex = lookupCstPool(atomSigPool, eachAtomSigs[i]);
00679         if(poolIndex==-1)
00680         {
00681             atomSigPool.push_back(eachAtomSigs[i]);
00682             poolIndex = atomSigPool.size()-1;
00683         }
00684         atomData[i].atomSigIdx = poolIndex;
00685     }
00686 
00687     printf("Atom's sigs: %d\n", (int)atomSigPool.size());
00688 
00689     delete[] eachAtomSigs;
00690 }
00691 
00692 void outputPsfFile(FILE *ofp)
00693 {
00694 
00695     fprintf(ofp, "%d !NSEGMENTNAMES\n", segNamePool.size());
00696     for(int i=0; i<segNamePool.size(); i++)
00697     {
00698         fprintf(ofp, "%s\n", segNamePool[i].c_str());
00699     }
00700 
00701     fprintf(ofp, "%d !NRESIDUENAMES\n", resNamePool.size());
00702     for(int i=0; i<resNamePool.size(); i++)
00703     {
00704         fprintf(ofp, "%s\n", resNamePool[i].c_str());
00705     }
00706 
00707     fprintf(ofp, "%d !NATOMNAMES\n", atomNamePool.size());
00708     for(int i=0; i<atomNamePool.size(); i++)
00709     {
00710         fprintf(ofp, "%s\n", atomNamePool[i].c_str());
00711     }
00712 
00713     fprintf(ofp, "%d !NATOMTYPES\n", atomTypePool.size());
00714     for(int i=0; i<atomTypePool.size(); i++)
00715     {
00716         fprintf(ofp, "%s\n", atomTypePool[i].c_str());
00717     }
00718 
00719     fprintf(ofp, "%d !NCHARGES\n", chargePool.size());
00720     for(int i=0; i<chargePool.size(); i++)
00721     {
00722         fprintf(ofp, "%f\n", chargePool[i]);
00723     }
00724 
00725     fprintf(ofp, "%d !NMASSES\n", massPool.size());
00726     for(int i=0; i<massPool.size(); i++)
00727     {
00728         fprintf(ofp, "%f\n", massPool[i]);
00729     }
00730 
00731 
00732     fprintf(ofp, "%d !NATOMSIGS\n", atomSigPool.size());
00733     for(int i=0; i<atomSigPool.size(); i++)
00734     {
00735         AtomSigInfo& oneAtomSig = atomSigPool[i];
00736         int oneTypeCnt = oneAtomSig.bondSigIndices.size();
00737         fprintf(ofp, "%d !%sSIGS\n", oneTypeCnt, "NBOND");
00738         for(int j=0; j<oneTypeCnt; j++)
00739         {
00740             short idx = oneAtomSig.bondSigIndices[j];
00741             TupleSignature& tSig = sigsOfBonds[idx];
00742             tSig.output(ofp);
00743         }
00744 
00745         oneTypeCnt = oneAtomSig.angleSigIndices.size();
00746         fprintf(ofp, "%d !%sSIGS\n", oneTypeCnt, "NTHETA");
00747         for(int j=0; j<oneTypeCnt; j++)
00748         {
00749             short idx = oneAtomSig.angleSigIndices[j];
00750             TupleSignature& tSig = sigsOfAngles[idx];
00751             tSig.output(ofp);
00752         }
00753 
00754         oneTypeCnt = oneAtomSig.dihedralSigIndices.size();
00755         fprintf(ofp, "%d !%sSIGS\n", oneTypeCnt, "NPHI");
00756         for(int j=0; j<oneTypeCnt; j++)
00757         {
00758             short idx = oneAtomSig.dihedralSigIndices[j];
00759             TupleSignature& tSig = sigsOfDihedrals[idx];
00760             tSig.output(ofp);
00761         }
00762 
00763         oneTypeCnt = oneAtomSig.improperSigIndices.size();
00764         fprintf(ofp, "%d !%sSIGS\n", oneTypeCnt, "NIMPHI");
00765         for(int j=0; j<oneTypeCnt; j++)
00766         {
00767             short idx = oneAtomSig.improperSigIndices[j];
00768             TupleSignature& tSig = sigsOfImpropers[idx];
00769             tSig.output(ofp);
00770         }
00771 
00772         oneTypeCnt = oneAtomSig.crosstermSigIndices.size();
00773         fprintf(ofp, "%d !%sSIGS\n", oneTypeCnt, "NCRTERM");
00774         for(int j=0; j<oneTypeCnt; j++)
00775         {
00776             short idx = oneAtomSig.crosstermSigIndices[j];
00777             TupleSignature& tSig = sigsOfCrossterms[idx];
00778             tSig.output(ofp);
00779         }
00780     }
00781 
00782     //2. Output exclusion signatures
00783     int exclSigCnt = sigsOfExclusions.size();
00784     fprintf(ofp, "%d !NEXCLSIGS\n", exclSigCnt);
00785     for(int i=0; i<exclSigCnt; i++)
00786     {
00787         ExclSigInfo *sig = &sigsOfExclusions[i];
00788         //first line is for full exclusions (1-2, 1-3) in the format of count offset1 offset2 offset3 ...
00789         fprintf(ofp, "%d", sig->fullExclOffset.size());
00790         for(int j=0; j<sig->fullExclOffset.size(); j++)
00791             fprintf(ofp, " %d", sig->fullExclOffset[j]);
00792         fprintf(ofp, "\n");
00793 
00794         //second line is for modified exclusions (1-4)
00795         fprintf(ofp, "%d", sig->modExclOffset.size());
00796         for(int j=0; j<sig->modExclOffset.size(); j++)
00797             fprintf(ofp, " %d", sig->modExclOffset[j]);
00798         fprintf(ofp, "\n");
00799     }
00800 
00801     //3. Output atom info
00802     fprintf(ofp, "%d !NATOM\n", g_mol->numAtoms);
00803     for(int i=0; i<g_mol->numAtoms; i++)
00804     {
00805         BasicAtomInfo &one = atomData[i];
00806         fprintf(ofp, "%d %d %d %d %d %d %d %d %d %d\n", one.segNameIdx, one.resID, one.resNameIdx, one.atomNameIdx,
00807                 one.atomTypeIdx, one.chargeIdx, one.massIdx, one.atomSigIdx, one.exclSigIdx, eachAtomClusterID[i]);
00808     }
00809     //fprintf(ofp, "\n");
00810 
00811     delete[] atomData;
00812     delete[] eachAtomClusterID;
00813 
00814     //4.Output the parameter new values if extraBonds are present.
00815     if(g_simParam->extraBondsOn)
00816     {
00817         // 1) output extra params for bonds
00818         int numExtraParams = extraBondParams.size();
00819         fprintf(ofp, "%d!NEXTRABONDPARAMS\n", numExtraParams);
00820         if(numExtraParams>0)
00821         {
00822             vector<BondValue>::iterator bpIter;
00823             for(bpIter=extraBondParams.begin(); bpIter!=extraBondParams.end(); bpIter++)
00824             {
00825                 fprintf(ofp, "%f %f\n", bpIter->k, bpIter->x0);
00826             }
00827         }
00828 
00829         // 2) output extra params for angles
00830         numExtraParams = extraAngleParams.size();
00831         fprintf(ofp, "%d!NEXTRAANGLEPARAMS\n", numExtraParams);
00832         if(numExtraParams>0)
00833         {
00834             NAMD_die("Output extra angle params not implemented!");
00835         }
00836 
00837         // 3) output extra params for dihedrals
00838         numExtraParams = extraDihedralParams.size();
00839         fprintf(ofp, "%d!NEXTRADIHEDRALPARAMS\n", numExtraParams);
00840         if(numExtraParams>0)
00841         {
00842             NAMD_die("Output extra dihedral params not implemented!");
00843         }
00844 
00845         // 4) output extra params for impropers
00846         numExtraParams = extraImproperParams.size();
00847         fprintf(ofp, "%d!NEXTRAIMPROPERPARAMS\n", numExtraParams);
00848         if(numExtraParams>0)
00849         {
00850             NAMD_die("Output extra improper params not implemented!");
00851         }
00852     }
00853 
00854 
00855     //5. Output the "multiplicity" field TUPLE_array of the Parameter object
00856     fprintf(ofp, "!DIHEDRALPARAMARRAY\n");
00857     for(int i=0; i<g_param->NumDihedralParams; i++)
00858     {
00859         fprintf(ofp, "%d ", g_param->dihedral_array[i].multiplicity);
00860     }
00861     fprintf(ofp, "\n");
00862     fprintf(ofp, "!IMPROPERPARAMARRAY\n");
00863     for(int i=0; i<g_param->NumImproperParams; i++)
00864     {
00865         fprintf(ofp, "%d ", g_param->improper_array[i].multiplicity);
00866     }
00867     fprintf(ofp, "\n");
00868 }
00869 
00870 void getAtomData(FILE *ifp)
00871 {
00872     char buffer[512]; //Buffer for reading from file
00873     int atomID=0;
00874     int lastAtomID=0; //to ensure atoms are in order
00875     char segmentName[11];
00876     int residueID;
00877     char residueName[11];
00878     char atomName[11];
00879     char atomType[11];
00880     Real charge;
00881     Real mass;
00882     int fieldsCnt; //Number of fields read by sscanf
00883 
00884     int numAtoms = g_mol->numAtoms;
00885 
00886     //1. parse atom data to build constant pool (atom name, mass, charge etc.)
00887     atomData = new BasicAtomInfo[numAtoms];
00888     while(atomID < numAtoms)
00889     {
00890         NAMD_read_line(ifp, buffer);
00891         //If it is blank or a comment, skip it
00892         if(NAMD_blank_string(buffer) || buffer[0]=='!')
00893             continue;
00894 
00895         //Fields are arranged as follows:
00896         //atomID, segName, resID, resName, atomName, atomType, charge, mass
00897         fieldsCnt = sscanf(buffer, "%d %s %d %s %s %s %f %f",
00898                            &atomID, segmentName, &residueID, residueName,
00899                            atomName, atomType, &charge, &mass);
00900 
00901         //check to make sure we found what we were expecting
00902         if(fieldsCnt!=8)
00903         {
00904             char errMsg[128];
00905             sprintf(errMsg, "BAD ATOM LINE FORMAT IN PSF FILE FOR ATOM %d (%s)", lastAtomID+1, buffer);
00906             NAMD_die(errMsg);
00907         }
00908 
00909         // check if this is in XPLOR format
00910         int atomTypeNum;
00911         if(sscanf(atomType, "%d", &atomTypeNum)>0)
00912             NAMD_die("Structure (psf) file is either in CHARMM format (with numbers for atoms types, the X-PLOR format using names is required) or the segment name field is empty.");
00913 
00914         //make sure the atoms were in sequence
00915         if(atomID!=lastAtomID+1)
00916         {
00917             char errMsg[128];
00918             sprintf(errMsg, "ATOM NUMBERS OUT OF ORDER AT ATOM #%d IN FSF FILE",
00919                     lastAtomID+1);
00920             NAMD_die(errMsg);
00921         }
00922 
00923         lastAtomID++;
00924 
00925         //building constant pool
00926         int poolIndex;
00927         string fieldName;
00928         fieldName.assign(segmentName);
00929         poolIndex = lookupCstPool(segNamePool, fieldName);
00930         if(poolIndex==-1)
00931         {
00932             segNamePool.push_back(fieldName);
00933             poolIndex = segNamePool.size()-1;
00934         }
00935         atomData[atomID-1].segNameIdx = poolIndex;
00936 
00937         //poolIndex = lookupCstPool(resIDPool, residueID);
00938         //if(poolIndex==-1)
00939         //    resIDPool.push_back(residueID);
00940         atomData[atomID-1].resID = residueID;
00941 
00942         fieldName.assign(residueName);
00943         poolIndex = lookupCstPool(resNamePool, fieldName);
00944         if(poolIndex==-1)
00945         {
00946             resNamePool.push_back(fieldName);
00947             poolIndex = resNamePool.size()-1;
00948         }
00949         atomData[atomID-1].resNameIdx = poolIndex;
00950 
00951         fieldName.assign(atomName);
00952         poolIndex = lookupCstPool(atomNamePool, fieldName);
00953         if(poolIndex==-1)
00954         {
00955             atomNamePool.push_back(fieldName);
00956             poolIndex = atomNamePool.size()-1;
00957         }
00958         atomData[atomID-1].atomNameIdx = poolIndex;
00959 
00960         fieldName.assign(atomType);
00961         poolIndex = lookupCstPool(atomTypePool, fieldName);
00962         if(poolIndex==-1)
00963         {
00964             atomTypePool.push_back(fieldName);
00965             poolIndex = atomTypePool.size()-1;
00966         }
00967         atomData[atomID-1].atomTypeIdx = poolIndex;
00968 
00969         poolIndex = lookupCstPool(chargePool, charge);
00970         if(poolIndex==-1)
00971         {
00972             chargePool.push_back(charge);
00973             poolIndex = chargePool.size()-1;
00974         }
00975         atomData[atomID-1].chargeIdx = poolIndex;
00976 
00977         poolIndex = lookupCstPool(massPool, mass);
00978         if(poolIndex==-1)
00979         {
00980             massPool.push_back(mass);
00981             poolIndex = massPool.size()-1;
00982         }
00983         atomData[atomID-1].massIdx = poolIndex;
00984     }
00985 }
00986 
00987 void getExtraBonds(StringList *file)
00988 {
00989     char err_msg[512];
00990     int a1, a2, a3, a4;
00991     float k, ref;
00992 
00993     int numAtoms = g_mol->numAtoms;
00994 
00995     if(!file)
00996     {
00997         NAMD_die("NO EXTRA BONDS FILES SPECIFIED");
00998     }
00999 
01000     for ( ; file; file = file->next )
01001     {  // loop over files
01002         FILE *f = fopen(file->data,"r");
01003         if ( ! f )
01004         {
01005             sprintf(err_msg, "UNABLE TO OPEN EXTRA BONDS FILE %s", file->data);
01006             NAMD_err(err_msg);
01007         }
01008         else
01009         {
01010             iout << iINFO << "READING EXTRA BONDS FILE " << file->data <<"\n"<<endi;
01011         }
01012 
01013         while ( 1 )
01014         {
01015             char buffer[512];
01016             int ret_code;
01017             int poolIndex;
01018             do
01019             {
01020                 ret_code = NAMD_read_line(f, buffer);
01021             }
01022             while ( (ret_code==0) && (NAMD_blank_string(buffer)) );
01023             if (ret_code!=0)
01024                 break;
01025 
01026             char type[512];
01027             sscanf(buffer,"%s",type);
01028 
01029 #define CHECKATOMID(ATOMID) if ( ATOMID < 0 || ATOMID >= numAtoms ) badatom = 1;
01030 
01031             int badline = 0;
01032             int badatom = 0;
01033             if ( ! strncasecmp(type,"bond",4) )
01034             {
01035                 if ( sscanf(buffer, "%s %d %d %f %f %s",
01036                             type, &a1, &a2, &k, &ref, err_msg) != 5 )
01037                     badline = 1;
01038                 else
01039                 {
01040                     CHECKATOMID(a1)
01041                     CHECKATOMID(a2)
01042                 }
01043                 Bond tmp;                
01044                 tmp.atom1 = a1;
01045                 tmp.atom2 = a2;
01046 
01047                 BondValue tmpv;
01048                 tmpv.k = k;
01049                 tmpv.x0 = ref;
01050 
01051                 poolIndex = lookupCstPool(extraBondParams, tmpv);
01052                 if(poolIndex==-1){               
01053                     extraBondParams.push_back(tmpv);
01054                     poolIndex = extraBondParams.size()-1;
01055                 }                 
01056                 tmp.bond_type = g_param->NumBondParams + poolIndex;
01057                 extraBonds.push_back(tmp);
01058             }
01059             else if ( ! strncasecmp(type,"angle",4) )
01060             {
01061                 if ( sscanf(buffer, "%s %d %d %d %f %f %s",
01062                             type, &a1, &a2, &a3, &k, &ref, err_msg) != 6 )
01063                     badline = 1;
01064                 else
01065                 {
01066                     CHECKATOMID(a1)
01067                     CHECKATOMID(a2)
01068                     CHECKATOMID(a3)
01069                 }
01070                 Angle tmp;
01071                 tmp.atom1 = a1;
01072                 tmp.atom2 = a2;
01073                 tmp.atom3 = a3;
01074                 
01075                 AngleValue tmpv;
01076                 tmpv.k = k;
01077                 tmpv.theta0 = ref / 180. * PI;
01078                 tmpv.k_ub = 0;
01079                 tmpv.r_ub = 0;
01080 
01081                 poolIndex = lookupCstPool(extraAngleParams, tmpv);
01082                 if(poolIndex==-1){               
01083                     extraAngleParams.push_back(tmpv);
01084                     poolIndex = extraAngleParams.size()-1;
01085                 }                 
01086                 tmp.angle_type = g_param->NumAngleParams + poolIndex;
01087                 extraAngles.push_back(tmp);
01088             }
01089             else if ( ! strncasecmp(type,"dihedral",4) )
01090             {
01091                 if ( sscanf(buffer, "%s %d %d %d %d %f %f %s",
01092                             type, &a1, &a2, &a3, &a4, &k, &ref, err_msg) != 7 )
01093                     badline = 1;
01094                 else
01095                 {
01096                     CHECKATOMID(a1)
01097                     CHECKATOMID(a2)
01098                     CHECKATOMID(a3)
01099                     CHECKATOMID(a4)
01100                 }
01101                 Dihedral tmp;
01102                 tmp.atom1 = a1;
01103                 tmp.atom2 = a2;
01104                 tmp.atom3 = a3;
01105                 tmp.atom4 = a4;
01106                 
01107                 DihedralValue tmpv;
01108                 tmpv.multiplicity = 1;
01109                 tmpv.values[0].n = 0;
01110                 tmpv.values[0].k = k;
01111                 tmpv.values[0].delta = ref / 180. * PI;
01112 
01113                 poolIndex = lookupCstPool(extraDihedralParams, tmpv);
01114                 if(poolIndex==-1){               
01115                     extraDihedralParams.push_back(tmpv);
01116                     poolIndex = extraDihedralParams.size()-1;
01117                 }                 
01118                 tmp.dihedral_type = g_param->NumDihedralParams + poolIndex;
01119                 extraDihedrals.push_back(tmp);
01120             }
01121             else if ( ! strncasecmp(type,"improper",4) )
01122             {
01123                 if ( sscanf(buffer, "%s %d %d %d %d %f %f %s",
01124                             type, &a1, &a2, &a3, &a4, &k, &ref, err_msg) != 7 )
01125                     badline = 1;
01126                 else
01127                 {
01128                     CHECKATOMID(a1)
01129                     CHECKATOMID(a2)
01130                     CHECKATOMID(a3)
01131                     CHECKATOMID(a4)
01132                 }
01133                 Improper tmp;
01134                 tmp.atom1 = a1;
01135                 tmp.atom2 = a2;
01136                 tmp.atom3 = a3;
01137                 tmp.atom4 = a4;
01138                 
01139                 ImproperValue tmpv;
01140                 tmpv.multiplicity = 1;
01141                 tmpv.values[0].n = 0;
01142                 tmpv.values[0].k = k;
01143                 tmpv.values[0].delta = ref / 180. * PI;
01144 
01145                 poolIndex = lookupCstPool(extraImproperParams, tmpv);
01146                 if(poolIndex==-1){               
01147                     extraImproperParams.push_back(tmpv);
01148                     poolIndex = extraImproperParams.size()-1;
01149                 }                 
01150                 tmp.improper_type = g_param->NumImproperParams + poolIndex;
01151                 extraImpropers.push_back(tmp);
01152             }
01153             else if ( ! strncasecmp(type,"#",1) )
01154             {
01155                 continue;  // comment
01156             }
01157             else
01158             {
01159                 badline = 1;
01160             }
01161 #undef CHECKATOMID
01162             if ( badline )
01163             {
01164                 sprintf(err_msg, "BAD LINE IN EXTRA BONDS FILE %s: %s",
01165                         file->data, buffer);
01166                 NAMD_die(err_msg);
01167             }
01168             if ( badatom )
01169             {
01170                 sprintf(err_msg, "BAD ATOM ID IN EXTRA BONDS FILE %s: %s",
01171                         file->data, buffer);
01172                 NAMD_die(err_msg);
01173             }
01174         }
01175         fclose(f);
01176     }  // loop over files
01177 }
01178 
01179 
01180 //In this function, safety check will also be performed:
01181 //1. bond itself 2.duplicate bonds such as (a, b) & (b, a)
01182 void getBondData(FILE *fd)
01183 {
01184 
01185     //first read bonds
01186     int atom_nums[2];  // Atom indexes for the bonded atoms
01187     char atom1name[11];  // Atom type for atom #1
01188     char atom2name[11];  // Atom type for atom #2
01189     register int j;      // Loop counter
01190     int num_read=0;    // Number of bonds read so far
01191 
01192     int numBonds = g_mol->numBonds;
01193     int origNumBonds = numBonds;   // number of bonds in file header
01194 
01195     /*  Allocate the array to hold the bonds      */
01196     Bond *bonds=new Bond[numBonds + extraBonds.size()];
01197 
01198     if (bonds == NULL)
01199     {
01200         NAMD_die("memory allocations failed in Molecule::read_bonds");
01201     }
01202 
01203     /*  Loop through and read in all the bonds      */
01204     while (num_read < numBonds)
01205     {
01206         /*  Loop and read in the two atom indexes    */
01207         for (j=0; j<2; j++)
01208         {
01209             /*  Read the atom number from the file.         */
01210             /*  Subtract 1 to convert the index from the    */
01211             /*  1 to NumAtoms used in the file to the       */
01212             /*  0 to NumAtoms-1 that we need    */
01213             atom_nums[j]=NAMD_read_int(fd, "BONDS")-1;
01214 
01215             /*  Check to make sure the index isn't too big  */
01216             if (atom_nums[j] >= g_mol->numAtoms)
01217             {
01218                 char err_msg[128];
01219 
01220                 sprintf(err_msg, "BOND INDEX %d GREATER THAN NATOM %d IN BOND # %d IN PSF FILE", atom_nums[j]+1, g_mol->numAtoms, num_read+1);
01221                 NAMD_die(err_msg);
01222             }
01223         }
01224 
01225         if(atom_nums[0] == atom_nums[1])
01226         { //bond to itself
01227             char err_msg[128];
01228             sprintf(err_msg, "ATOM %d is bonded to itself!", atom_nums[0]+1);
01229             NAMD_die(err_msg);
01230         }
01231 
01232         /*  Get the atom type for the two atoms.  When we query */
01233         /*  the parameter object, we need to send the atom type */
01234         /*  that is alphabetically first as atom 1.    */
01235         const char *atom1Type = atomTypePool[atomData[atom_nums[0]].atomTypeIdx].c_str();
01236         const char *atom2Type = atomTypePool[atomData[atom_nums[1]].atomTypeIdx].c_str();
01237         if (strcasecmp(atom1Type,atom2Type) < 0)
01238         {
01239             strcpy(atom1name, atom1Type);
01240             strcpy(atom2name, atom2Type);
01241         }
01242         else
01243         {
01244             strcpy(atom2name, atom1Type);
01245             strcpy(atom1name, atom2Type);
01246         }
01247 
01248         /*  Query the parameter object for the constants for    */
01249         /*  this bond            */
01250         Bond *b = &(bonds[num_read]);
01251         g_param->assign_bond_index(atom1name, atom2name, b);
01252 
01253         /*  Assign the atom indexes to the array element  */
01254         b->atom1=atom_nums[0];
01255         b->atom2=atom_nums[1];
01256 
01257         /*  Make sure this isn't a fake bond meant for shake in x-plor.  */
01258         Real k, x0;
01259         g_param->get_bond_params(&k,&x0,b->bond_type);
01260         if ( k == 0. )
01261             --numBonds;  // fake bond
01262         else
01263             ++num_read;  // real bond
01264     }
01265 
01266     /*  Tell user about our subterfuge  */
01267     if ( numBonds != origNumBonds )
01268     {
01269         iout << iWARN << "Ignored " << origNumBonds - numBonds <<
01270         " bonds with zero force constants.\n" << endi;
01271         iout << iWARN <<
01272         "Will get H-H distance in rigid H2O from H-O-H angle.\n" << endi;
01273     }
01274 
01275     //copy extra bonds to bonds structure
01276     int numRealBonds = numBonds;
01277     for(int i=0; i<extraBonds.size(); i++)
01278         bonds[numBonds+i] = extraBonds[i];
01279     numBonds += extraBonds.size();
01280     extraBonds.clear();
01281 
01282     g_mol->numBonds = numBonds;
01283 
01284     //then creating bond's tupleSignature
01285     for(int i=0; i<numBonds; i++)
01286     {
01287         Bond *b = bonds+i;
01288         TupleSignature oneSig(1,BOND,b->bond_type);
01289         oneSig.offset[0] = b->atom2 - b->atom1;
01290         oneSig.isReal = (i<numRealBonds );
01291 
01292         int poolIndex = lookupCstPool(sigsOfBonds, oneSig);
01293         int newSig=0;
01294         if(poolIndex == -1)
01295         {
01296             sigsOfBonds.push_back(oneSig);
01297             poolIndex = (short)sigsOfBonds.size()-1;
01298             newSig=1;
01299         }
01300 
01301         if(!newSig)
01302         {//check duplicate bonds in the form of (a, b) && (a, b);
01303             int dupIdx = lookupCstPool(eachAtomSigs[b->atom1].bondSigIndices, (short)poolIndex);
01304             if(dupIdx!=-1)
01305             {
01306                 char err_msg[128];
01307                 sprintf(err_msg, "Duplicate bond %d-%d!", b->atom1+1, b->atom2+1);
01308                 NAMD_die(err_msg);
01309             }
01310         }
01311         eachAtomSigs[b->atom1].bondSigIndices.push_back(poolIndex);
01312     }
01313 
01314     //check duplicate bonds in the form of (a, b) && (b, a)
01315     for(int i=0; i<numBonds; i++)
01316     {
01317         Bond *b=bonds+i;
01318         int atom2 = b->atom2;
01319         int thisOffset = atom2 - b->atom1;
01320         for(int j=0; j<eachAtomSigs[atom2].bondSigIndices.size(); j++)
01321         {
01322             short atom2BondId = eachAtomSigs[atom2].bondSigIndices[j];
01323             TupleSignature *secSig = &(sigsOfBonds[atom2BondId]);
01324             if(thisOffset== -(secSig->offset[0]))
01325             {
01326                 char err_msg[128];
01327                 sprintf(err_msg, "Duplicate bond %d-%d because two atoms are just reversed!", b->atom1+1, atom2+1);
01328                 NAMD_die(err_msg);
01329             }
01330         }
01331     }      
01332 
01333     //building clusters for this simulation system in two steps
01334     //1. create a list for each atom where each atom in the list is bonded with that atom
01335     vector<int> *atomListOfBonded = new vector<int>[g_mol->numAtoms];
01336 
01337     for(int i=0; i<numRealBonds; i++)
01338     {
01339         Bond *b=bonds+i;
01340         int atom1 = b->atom1;
01341         int atom2 = b->atom2;
01342         atomListOfBonded[atom1].push_back(atom2);
01343         atomListOfBonded[atom2].push_back(atom1);
01344     }
01345 
01346     delete [] bonds;
01347 
01348     //2. using breadth-first-search to build the clusters. Here, we avoid recursive call
01349     // because the depth of calls may be of thousands which will blow up the stack, and
01350     //recursive call is slower than the stack-based BFS.
01351     //Considering such structure
01352     //1->1245; 7->1243; 1243->1245
01353     eachAtomClusterID = new int[g_mol->numAtoms];
01354     for(int i=0; i<g_mol->numAtoms; i++)
01355         eachAtomClusterID[i] = -1;
01356 
01357     for(int i=0; i<g_mol->numAtoms; i++)
01358     {
01359         int curClusterID=eachAtomClusterID[i];
01360         if(curClusterID==-1)
01361         {
01362             curClusterID=i;
01363         }
01364 
01365         deque<int> toVisitAtoms;
01366         toVisitAtoms.push_back(i);
01367         while(!toVisitAtoms.empty())
01368         {
01369             int visAtomID = toVisitAtoms.front();
01370             toVisitAtoms.pop_front();
01371             eachAtomClusterID[visAtomID] = curClusterID;
01372             for(int j=0; j<atomListOfBonded[visAtomID].size(); j++)
01373             {
01374                 int otherAtom = atomListOfBonded[visAtomID][j];
01375                 if(eachAtomClusterID[otherAtom]!=curClusterID)
01376                     toVisitAtoms.push_back(otherAtom);
01377             }
01378         }
01379     }
01380 
01381     //Now the clusterID of each atom should be in the non-decreasing order, otherwise
01382     //report a psf format error. Well, it's not a real problem, it can be fixed the problem
01383     //later. Here, such assumption is made for the ease of programming
01384     int curClusterID;
01385     int prevClusterID=eachAtomClusterID[0];
01386     int curClusterSize=1;
01387     for(int i=1; i<g_mol->numAtoms; i++)
01388     {
01389         curClusterID = eachAtomClusterID[i];
01390         if(curClusterID > prevClusterID)
01391         {
01392             eachClusterSize.push_back(curClusterSize);
01393             curClusterSize=1;
01394         }
01395         else if(curClusterID == prevClusterID)
01396         {
01397             curClusterSize++;
01398         }
01399         else
01400         { //psf format error
01401             char errMsg[256];
01402             sprintf(errMsg, "Cluster ID of each atom should be in increasing order (error for atom %d)!\n", i+1);
01403             NAMD_die(errMsg);
01404         }
01405         prevClusterID = curClusterID;
01406     }
01407     eachClusterSize.push_back(curClusterSize); //record the last sequence of cluster
01408 
01409     //Now iterate over the cluster size again to filter out the repeating cluster size.
01410     //Since the cluster id of atoms is monotonically increasing, the size of the cluster can be used
01411     //as this cluster's signature.
01412     //After this, eachAtomClusterID will store the cluster signature. Only the atom whose id is "clustersize"
01413     //will store the size. Others will be -1
01414 
01415     int aid=0;
01416     for(int clusterIdx=0; clusterIdx<eachClusterSize.size(); clusterIdx++)
01417     {
01418         int curSize = eachClusterSize[clusterIdx];
01419         eachAtomClusterID[aid] = curSize;
01420         for(int i=aid+1; i<aid+curSize; i++)
01421             eachAtomClusterID[i] = -1;
01422         aid += curSize;
01423     }
01424 
01425     //check whether cluster is built correctly
01426     /*printf("num clusters: %d\n", eachClusterSize.size());
01427     FILE *checkFile = fopen("cluster.opt", "w");
01428     for(int i=0; i<g_mol->numAtoms; i++)  fprintf(checkFile, "%d\n", eachAtomClusterID[i]);
01429     fclose(checkFile);*/
01430 
01431     eachClusterSize.clear();
01432     for(int i=0; i<g_mol->numAtoms; i++)
01433         atomListOfBonded[i].clear();
01434     delete [] atomListOfBonded;
01435 }
01436 
01437 void getAngleData(FILE *fd)
01438 {
01439 
01440     int atom_nums[3];  //  Atom numbers for the three atoms
01441     char atom1name[11];  //  Atom type for atom 1
01442     char atom2name[11];  //  Atom type for atom 2
01443     char atom3name[11];  //  Atom type for atom 3
01444     register int j;      //  Loop counter
01445     int num_read=0;    //  Number of angles read so far
01446 
01447     int numAngles = g_mol->numAngles;
01448     int origNumAngles = numAngles;  // Number of angles in file
01449     /*  Alloc the array of angles          */
01450     Angle *angles=new Angle[numAngles + extraAngles.size()];
01451 
01452     if (angles == NULL)
01453     {
01454         NAMD_die("memory allocation failed in Molecule::read_angles");
01455     }
01456 
01457     /*  Loop through and read all the angles      */
01458     while (num_read < numAngles)
01459     {
01460         /*  Loop through the 3 atom indexes in the current angle*/
01461         for (j=0; j<3; j++)
01462         {
01463             /*  Read the atom number from the file.         */
01464             /*  Subtract 1 to convert the index from the    */
01465             /*  1 to NumAtoms used in the file to the       */
01466             /*  0 to NumAtoms-1 that we need    */
01467             atom_nums[j]=NAMD_read_int(fd, "ANGLES")-1;
01468 
01469             /*  Check to make sure the atom index doesn't   */
01470             /*  exceed the Number of Atoms      */
01471             if (atom_nums[j] >= g_mol->numAtoms)
01472             {
01473                 char err_msg[128];
01474 
01475                 sprintf(err_msg, "ANGLES INDEX %d GREATER THAN NATOM %d IN ANGLES # %d IN PSF FILE", atom_nums[j]+1, g_mol->numAtoms, num_read+1);
01476                 NAMD_die(err_msg);
01477             }
01478         }
01479 
01480         /*  Place the bond name that is alphabetically first  */
01481         /*  in the atom1name.  This is OK since the order of    */
01482         /*  atom1 and atom3 are interchangable.  And to search  */
01483         /*  the tree of angle parameters, we need the order     */
01484         /*  to be predictable.          */
01485 
01486         const char *atom1Type = atomTypePool[atomData[atom_nums[0]].atomTypeIdx].c_str();
01487         const char *atom2Type = atomTypePool[atomData[atom_nums[1]].atomTypeIdx].c_str();
01488         const char *atom3Type = atomTypePool[atomData[atom_nums[2]].atomTypeIdx].c_str();
01489 
01490         if (strcasecmp(atom1Type,atom2Type) < 0)
01491         {
01492             strcpy(atom1name, atom1Type);
01493             strcpy(atom2name, atom2Type);
01494             strcpy(atom3name, atom3Type);
01495         }
01496         else
01497         {
01498             strcpy(atom1name, atom1Type);
01499             strcpy(atom2name, atom2Type);
01500             strcpy(atom3name, atom3Type);
01501         }
01502 
01503         /*  Get the constant values for this bond from the  */
01504         /*  parameter object          */
01505         g_param->assign_angle_index(atom1name, atom2name,
01506                                     atom3name, &(angles[num_read]));
01507 
01508         /*  Assign the three atom indices      */
01509         angles[num_read].atom1=atom_nums[0];
01510         angles[num_read].atom2=atom_nums[1];
01511         angles[num_read].atom3=atom_nums[2];
01512 
01513         /*  Make sure this isn't a fake angle meant for shake in x-plor.  */
01514         Real k, t0, k_ub, r_ub;
01515         g_param->get_angle_params(&k,&t0,&k_ub,&r_ub,angles[num_read].angle_type);
01516         if ( k == 0. && k_ub == 0. )
01517             --numAngles;  // fake angle
01518         else
01519             ++num_read;  // real angle
01520     }
01521 
01522     /*  Tell user about our subterfuge  */
01523     if ( numAngles != origNumAngles )
01524     {
01525         iout << iWARN << "Ignored " << origNumAngles - numAngles <<
01526         " angles with zero force constants.\n" << endi;
01527     }
01528 
01529     //copy extra angles to angles structure
01530     for(int i=0; i<extraAngles.size(); i++)
01531         angles[numAngles+i] = extraAngles[i];
01532     numAngles += extraAngles.size();
01533     extraAngles.clear();
01534 
01535     g_mol->numAngles = numAngles;
01536 
01537     //create angles' tupleSignature
01538     for(int i=0; i<numAngles; i++)
01539     {
01540         Angle *tuple = angles+i;
01541         TupleSignature oneSig(2,ANGLE,tuple->angle_type);
01542         int offset[2];
01543         offset[0] = tuple->atom2 - tuple->atom1;
01544         offset[1] = tuple->atom3 - tuple->atom1;
01545         oneSig.setOffsets(offset);
01546 
01547         int poolIndex = lookupCstPool(sigsOfAngles, oneSig);
01548         if(poolIndex == -1)
01549         {
01550             sigsOfAngles.push_back(oneSig);
01551             poolIndex = (short)sigsOfAngles.size()-1;
01552         }
01553         eachAtomSigs[tuple->atom1].angleSigIndices.push_back(poolIndex);
01554     }
01555 
01556     delete [] angles;
01557 
01558 }
01559 
01560 void getDihedralData(FILE *fd)
01561 {
01562 
01563     int atom_nums[4];  // The 4 atom indexes
01564     int last_atom_nums[4];  // Atom numbers from previous bond
01565     char atom1name[11];  // Atom type for atom 1
01566     char atom2name[11];  // Atom type for atom 2
01567     char atom3name[11];  // Atom type for atom 3
01568     char atom4name[11];  // Atom type for atom 4
01569     register int j;      // loop counter
01570     int num_read=0;    // number of dihedrals read so far
01571     int multiplicity=1;  // multiplicity of the current bond
01572     Bool duplicate_bond;  // Is this a duplicate of the last bond
01573     int num_unique=0;   // Number of unique dihedral bonds
01574 
01575     //  Initialize the array used to check for duplicate dihedrals
01576     for (j=0; j<4; j++)
01577         last_atom_nums[j] = -1;
01578 
01579     int numDihedrals = g_mol->numDihedrals;
01580     int numAtoms = g_mol->numAtoms;
01581 
01582     /*  Allocate an array to hold the Dihedrals      */
01583     Dihedral *dihedrals = new Dihedral[numDihedrals + extraDihedrals.size()];
01584 
01585     if (dihedrals == NULL)
01586     {
01587         NAMD_die("memory allocation failed in Molecule::read_dihedrals");
01588     }
01589 
01590     /*  Loop through and read all the dihedrals      */
01591     while (num_read < numDihedrals)
01592     {
01593         duplicate_bond = TRUE;
01594 
01595         /*  Loop through and read the 4 indexes for this bond   */
01596         for (j=0; j<4; j++)
01597         {
01598             /*  Read the atom number from the file.         */
01599             /*  Subtract 1 to convert the index from the    */
01600             /*  1 to NumAtoms used in the file to the       */
01601             /*  0 to NumAtoms-1 that we need    */
01602             atom_nums[j]=NAMD_read_int(fd, "DIHEDRALS")-1;
01603 
01604             /*  Check for an atom index that is too large  */
01605             if (atom_nums[j] >= numAtoms)
01606             {
01607                 char err_msg[128];
01608 
01609                 sprintf(err_msg, "DIHEDRALS INDEX %d GREATER THAN NATOM %d IN DIHEDRALS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
01610                 NAMD_die(err_msg);
01611             }
01612 
01613             //  Check to see if this atom matches the last bond
01614             if (atom_nums[j] != last_atom_nums[j])
016