PDB.C

Go to the documentation of this file.
00001 
00007 /*
00008    PDB Class
00009      Code to implement the PDB class.  This reads in a bunch of
00010    PDBData records from a file, given the filename.  See PDB.h
00011    for a bit more information.
00012 */
00013 
00014 #include <stdio.h>
00015 #include <stdlib.h>
00016 #ifndef WIN32
00017 #include <strings.h>
00018 #endif
00019 #include "common.h"
00020 #include "PDB.h"
00021 #include "SortableResizeArray.h"
00022 
00023 
00024 #ifdef MEM_OPT_VERSION
00025 //read in PDB from a file and eliminate the temporary memory usage of pdb atom list    
00026 PDB::PDB(const char *pdbfilename, int expectedNumAtoms){
00027   FILE *infile;
00028   char buf[160];
00029 
00030   altlocArray = 0;
00031   atomArray = new PDBCoreData[expectedNumAtoms];
00032 
00033   atomCount = 0;
00034   atomListHead = atomListTail = NULL;
00035   infile = Fopen(pdbfilename, "r");
00036   if (! infile) {
00037      char s[500];
00038      sprintf(s, "Cannot open file '%s' for input in PDB::PDB.", pdbfilename);
00039      NAMD_err(s);
00040   }
00041   
00042     // loop through each line in the file and get a record
00043   while ( fgets(buf, 150, infile) ) {
00044    PDBData *newelement;
00045    char *s;
00046    for (s=buf; *s && *s!='\n'; s++)  // chop off the '\n'
00047     ;
00048    *s = 0;
00049    if ( s == (buf + 149) ) {
00050      char s[500];
00051      sprintf( s, "Line too long in pdbfile %s:\n%s\n", pdbfilename, buf);
00052      NAMD_die(s);
00053    }
00054    *(s+1) = 0;  // just to be on the safe side
00055 
00056      // I now have a string; make a PDBData element out of it
00057    newelement = new_PDBData(buf);
00058    if (!newelement) {
00059       NAMD_die("Could not allocate PDBData.\n");
00060    }
00061      // I only know how to deal with ATOM and HETATM types; and
00062      //  I want to throw away the unknown data types; so
00063    if (newelement -> type() != PDBData::ATOM && 
00064            newelement -> type() != PDBData::HETATM) {
00065        delete newelement;
00066    } else {
00067        PDBAtom *thisAtom = (PDBAtom *)newelement;
00068        const BigReal *coor = thisAtom->coordinates();
00069        atomArray[atomCount].coor[0] = coor[0];
00070        atomArray[atomCount].coor[1] = coor[1];
00071        atomArray[atomCount].coor[2] = coor[2];
00072        atomArray[atomCount].myoccupancy = thisAtom->occupancy();
00073        atomArray[atomCount].tempfactor = thisAtom->temperaturefactor();
00074        atomCount++;
00075        delete newelement;
00076        if(atomCount > expectedNumAtoms) NAMD_die("Number of pdb and psf atoms are not the same!");       
00077    }
00078   }  // input while loop
00079  Fclose(infile);
00080 }
00081 #endif
00082 
00083 PDB::PDB(molfile_plugin_t *pIOHdl, void *pIOFileHdl, int numAtoms, const float *occupancy, const float *bfactor){
00084 #ifdef MEM_OPT_VERSION
00085   NAMD_die("Sorry, plugin IO is not supported in the memory optimized version.");
00086 #else
00087   molfile_timestep_t ts;
00088   float *atomcoords;
00089   memset(&ts, 0, sizeof(molfile_timestep_t));
00090 
00091   /* set defaults for unit cell information */
00092   ts.A = ts.B = ts.C = 0.0f;
00093   ts.alpha = ts.beta = ts.gamma = 90.0f; 
00094 
00095   atomCount = numAtoms;
00096 
00097   atomcoords = (float *) malloc(3*numAtoms*sizeof(float));
00098   memset(atomcoords, 0, 3*numAtoms*sizeof(float));
00099   ts.coords = atomcoords;
00100   
00101   if (pIOHdl->read_next_timestep(pIOFileHdl, numAtoms, &ts)) {
00102     free(atomcoords);    
00103     pIOHdl->close_file_read(pIOFileHdl);
00104     NAMD_die("ERROR: failed reading atom coordinates");
00105   }
00106 
00107   //load coordinates to PDB object
00108   //Note: the PDBAtom structure is very redundant in this case
00109   PDBAtom *tmpAtoms = new PDBAtom[numAtoms];
00110   atomArray = new PDBAtomPtr[numAtoms];
00111   BigReal tmpCoords[3];
00112   for(int i=0; i<numAtoms; i++) {
00113       PDBAtom *thisAtom = tmpAtoms+i;
00114       atomArray[i] = thisAtom;
00115 
00116       tmpCoords[0] = atomcoords[0];
00117       tmpCoords[1] = atomcoords[1];
00118       tmpCoords[2] = atomcoords[2];
00119       atomcoords += 3; //set to the next atom
00120       thisAtom->coordinates(tmpCoords);
00121   }
00122   free(ts.coords);
00123 
00124   if(occupancy) {
00125       for(int i=0; i<numAtoms; i++) {
00126           PDBAtom *thisAtom = tmpAtoms+i;
00127           thisAtom->occupancy((BigReal)occupancy[i]);          
00128       }
00129   }
00130   if(bfactor) {
00131       for(int i=0; i<numAtoms; i++) {
00132           PDBAtom *thisAtom = tmpAtoms+i;       
00133           thisAtom->temperaturefactor((BigReal)bfactor[i]);       
00134       }
00135   }
00136 #endif
00137 }
00138 
00139 
00140 // read in a file and stick all the elements on the appropriate list
00141 PDB::PDB( const char *pdbfilename) {
00142   FILE *infile;
00143   char buf[160];
00144 
00145   atomCount = 0;
00146   atomListHead = atomListTail = NULL;
00147   infile = Fopen(pdbfilename, "r");
00148   if (! infile) {
00149      char s[500];
00150      sprintf(s, "Cannot open file '%s' for input in PDB::PDB.", pdbfilename);
00151      NAMD_err(s);
00152   }
00153   
00154     // loop through each line in the file and get a record
00155   while ( fgets(buf, 150, infile) ) {
00156    PDBData *newelement;
00157    char *s;
00158    for (s=buf; *s && *s!='\n'; s++)  // chop off the '\n'
00159     ;
00160    *s = 0;
00161    if ( s == (buf + 149) ) {
00162      char s[500];
00163      sprintf( s, "Line too long in pdbfile %s:\n%s\n", pdbfilename, buf);
00164      NAMD_die(s);
00165    }
00166    *(s+1) = 0;  // just to be on the safe side
00167 
00168      // I now have a string; make a PDBData element out of it
00169    newelement = new_PDBData(buf);
00170    if (!newelement) {
00171       NAMD_die("Could not allocate PDBData.\n");
00172    }
00173      // I only know how to deal with ATOM and HETATM types; and
00174      //  I want to throw away the unknown data types; so
00175    if (newelement -> type() != PDBData::ATOM && 
00176            newelement -> type() != PDBData::HETATM) {
00177        delete newelement;
00178    } else {
00179        add_atom_element( (PDBAtom *) newelement);
00180    }
00181   }  // input while loop
00182  Fclose(infile);
00183  
00184  // now I have a linked list, and I know the size.  However,
00185  // that's got pretty slow access time for ramdom fetches, so
00186  // I'll turn it into an array
00187  {
00188 
00189 #ifdef MEM_OPT_VERSION
00190 //pruning the PDBData structure into PDBCoreData by only leaving X/Y/Z, occupancy and temperaturefactor
00191      altlocArray = new char[atomCount];
00192      atomArray = new PDBCoreData[atomCount];
00193      if(atomArray == NULL)
00194          NAMD_die("memory allocation failed in PDB::PDB");
00195 
00196      PDBAtomList *tmp = atomListHead;
00197      for(int i=0; tmp!=NULL; tmp=tmp->next, i++){
00198          const BigReal *coor = tmp->data->coordinates();
00199          atomArray[i].coor[0] = coor[0];
00200          atomArray[i].coor[1] = coor[1];
00201          atomArray[i].coor[2] = coor[2];
00202          atomArray[i].myoccupancy = tmp->data->occupancy();
00203          atomArray[i].tempfactor = tmp->data->temperaturefactor();
00204          altlocArray[i] = tmp->data->alternatelocation()[0];
00205      }
00206 
00207      //free the PDBAtomList
00208      tmp = atomListHead;
00209      while(tmp!=NULL){
00210          PDBAtomList *nextTmp = tmp->next;
00211          delete tmp->data;
00212          delete tmp;
00213          tmp = nextTmp;
00214      }
00215      atomListHead = atomListTail = NULL;
00216 #else
00217   atomArray = new PDBAtomPtr[atomCount];
00218   if ( atomArray == NULL )
00219   {
00220     NAMD_die("memory allocation failed in PDB::PDB");
00221   }
00222   PDBAtomList *tmp = atomListHead;
00223   int i=0;                              // just need to copy the pointers
00224   for (i=0, tmp = atomListHead; tmp != NULL; tmp = tmp -> next, i++) {
00225     atomArray[i] = tmp -> data;
00226   }
00227   // now delete the linked list (w/o deleting the data)
00228   PDBAtomList *tmp2;
00229   for (tmp2 = tmp = atomListHead; tmp != NULL; tmp = tmp2) {
00230     tmp2 = tmp->next;
00231     delete tmp;
00232   }
00233   atomListHead = atomListTail = NULL;
00234 #endif 
00235  }  // everything converted
00236  
00237 }
00238 
00239 //  Destructor - delete all the data pointed to by the array
00240 //   and then delete the array
00241 PDB::~PDB( void )
00242 {
00243 #ifndef MEM_OPT_VERSION
00244         int i;
00245         if ( atomArray[atomCount-1] == atomArray[0] + (atomCount-1) ) {
00246           delete [] atomArray[0];
00247         } else {
00248           for (i=atomCount-1; i>=0; i--)
00249             delete atomArray[i];
00250         }
00251 #else
00252         delete [] altlocArray;
00253 #endif
00254         delete [] atomArray;
00255         atomArray = NULL;
00256         atomCount = 0;
00257 }
00258 
00259 // print the PDB file out to a given file name
00260 void PDB::write(const char *outfilename, const char *commentline)
00261 {
00262         int i;
00263         char s[200];
00264         FILE *outfile;
00265         if ((outfile = fopen(outfilename, "w")) == NULL) {
00266            sprintf(s, "Cannot open file '%s' in PDB::write.", outfilename);
00267            NAMD_err(s);
00268         }
00269 
00270         if (commentline != NULL)
00271         {
00272                 sprintf(s, "REMARK  %s\n", commentline);
00273                 if (fputs(s, outfile) == EOF)
00274                 {
00275                         NAMD_err("EOF in PDB::write writing the comment line - file system full?");
00276                 }
00277         }
00278 
00279         for (i=0; i<atomCount; i++){ // I only contain ATOM/HETATM records
00280 #ifdef MEM_OPT_VERSION    
00281       atomArray[i].sprint(s, PDBData::COLUMNS);  
00282 #else
00283           atomArray[i]->sprint(s, PDBData::COLUMNS);
00284 #endif
00285           if ( (fputs(s, outfile)    == EOF) || 
00286                (fputs("\n", outfile) == EOF)    ) {
00287             sprintf(s, "EOF in PDB::write line %d - file system full?", i);
00288             NAMD_err(s);
00289           }
00290         }
00291         if (fputs("END\n", outfile) == EOF) {
00292            NAMD_err("EOF in PDB::write while printing 'END' -- file system full?");
00293         }
00294         if (fclose(outfile) == EOF) {
00295            NAMD_err("EOF in PDB::write while closing -- file system full?");
00296         }
00297           
00298 }
00299 
00300 // store the info on the linked list
00301 void PDB::add_atom_element( PDBAtom *newAtom)
00302 {
00303   PDBAtomList *tmp = new PDBAtomList;
00304   if ( tmp == NULL )
00305   {
00306     NAMD_die("memory allocation failed in PDB::add_atom_element");
00307   }
00308   tmp -> data = newAtom;
00309   tmp -> next = NULL;
00310   
00311   if (atomListHead == NULL) {        // make the list
00312     atomListHead = atomListTail = tmp;
00313   } else {
00314     atomListTail -> next = tmp;       // add to the tail
00315     atomListTail = tmp;
00316   }
00317   atomCount++;
00318 }
00319 
00320 
00321 // return the number of atoms found
00322 int PDB::num_atoms( void)
00323 {
00324   return atomCount;
00325 }
00326 
00327 
00328 // Reset all the atom positions.  This is used in preparation for
00329 // output in cases like the restart files, etc.
00330 void PDB::set_all_positions(Vector *pos)
00331 {
00332 #ifdef MEM_OPT_VERSION
00333     for(int i=0; i<atomCount; i++){
00334         atomArray[i].coor[0] = pos[i].x;
00335         atomArray[i].coor[1] = pos[i].y;
00336         atomArray[i].coor[2] = pos[i].z;
00337     }
00338 #else
00339         int i;
00340         PDBAtomPtr *atomptr;
00341 
00342         for (i=0, atomptr=atomArray; i<atomCount; atomptr++, i++)
00343         {
00344                 (*atomptr)->xcoor(pos[i].x);
00345                 (*atomptr)->ycoor(pos[i].y);
00346                 (*atomptr)->zcoor(pos[i].z);
00347         }
00348 #endif
00349 }
00350 
00351 void PDB::get_position_for_atom(Vector *pos, int aid){
00352 #ifdef MEM_OPT_VERSION
00353     pos->x = atomArray[aid].coor[0];
00354     pos->y = atomArray[aid].coor[1];
00355     pos->z = atomArray[aid].coor[2];
00356 #else
00357     PDBAtomPtr *atomptr = &atomArray[aid];
00358     pos->x = (*atomptr)->xcoor();
00359     pos->y = (*atomptr)->ycoor();
00360     pos->z = (*atomptr)->zcoor();
00361 #endif
00362 }
00363 
00364 //  Get all the atom positions into a list of Vectors
00365 void PDB::get_all_positions(Vector *pos)
00366 {
00367 #ifdef MEM_OPT_VERSION
00368     for(int i=0; i<atomCount; i++){
00369         pos[i].x = atomArray[i].coor[0];
00370         pos[i].y = atomArray[i].coor[1];
00371         pos[i].z = atomArray[i].coor[2];        
00372     }
00373 #else
00374         int i;
00375         PDBAtomPtr *atomptr;
00376 
00377         for (i=0, atomptr=atomArray; i<atomCount; atomptr++, i++)
00378         {
00379                 pos[i].x = (*atomptr)->xcoor();
00380                 pos[i].y = (*atomptr)->ycoor();
00381                 pos[i].z = (*atomptr)->zcoor();
00382         }
00383 #endif
00384 }
00385 
00386 //  given an index, return that atom
00387 #ifdef MEM_OPT_VERSION
00388 PDBCoreData *PDB::atom(int place){
00389     if(place<0 || place>=atomCount) return NULL;
00390     return &atomArray[place];
00391 }
00392 #else
00393 PDBAtom *PDB::atom(int place)
00394 {
00395   if (place <0 || place >= atomCount)
00396     return NULL;
00397   return atomArray[place];
00398 }
00399 #endif
00400 
00401 
00402 // find the lowest and highest bounds based on a fraction of the atoms
00403 //void PDB::find_extremes(BigReal *min, BigReal *max, Vector rec, BigReal frac) const
00404 void PDB::find_extremes_helper(
00405     SortableResizeArray<BigReal> &coor,
00406     BigReal &min, BigReal &max, Vector rec, BigReal frac
00407     )
00408 {
00409     SortableResizeArray<BigReal>::iterator c_i = coor.begin();
00410 #ifdef MEM_OPT_VERSION
00411     PDBCoreData *atomptr = atomArray;
00412     for(int i=0; i<atomCount; i++, atomptr++){
00413         Vector pos(atomptr->xcoor(),atomptr->ycoor(),atomptr->zcoor());
00414         c_i[i] = rec*pos;
00415     }
00416 #else
00417     PDBAtomPtr *atomptr = atomArray;
00418     for (int i=0; i<atomCount; ++i, ++atomptr) {
00419       PDBAtom *atom = *atomptr;
00420       Vector pos(atom->xcoor(),atom->ycoor(),atom->zcoor());
00421       c_i[i] = rec*pos;
00422     }
00423 #endif
00424     coor.sort();
00425     int ilow = (int)((1.0 - frac) * atomCount);
00426     if ( ilow < 0 ) ilow = 0;
00427     if ( ilow > atomCount/2 ) ilow = atomCount/2;
00428     int ihigh = atomCount - ilow - 1;
00429     BigReal span = coor[ihigh] - coor[ilow];
00430     BigReal extension = (1.0 - frac) * span / (2.0 * frac - 1.0);
00431     max = coor[ihigh] + extension;
00432     min = coor[ilow] - extension;
00433 }
00434 
00435 // Find the extreme edges of molecule in scaled coordinates,
00436 // where "frac" sets bounds based on a fraction of the atoms.
00437 void PDB::find_extremes(const Lattice &lattice, BigReal frac) {
00438   if (atomCount == 0) {
00439     smin = smax = 0;
00440   }
00441   else if (frac < 1.0) {
00442     // for now use the previous "sort the array" approach
00443     // for solving "frac"th largest and smallest selection problems
00444     SortableResizeArray<BigReal> coor;
00445     coor.resize(atomCount);  // allocate array space once
00446     find_extremes_helper(coor, smin.x, smax.x, lattice.a_r(), frac);
00447     find_extremes_helper(coor, smin.y, smax.y, lattice.b_r(), frac);
00448     find_extremes_helper(coor, smin.z, smax.z, lattice.c_r(), frac);
00449   }
00450   else {
00451     // finding absolute min and max does not require sorting
00452 #ifdef MEM_OPT_VERSION
00453     PDBCoreData *atomptr = atomArray;
00454     Vector p(atomptr->xcoor(),atomptr->ycoor(),atomptr->zcoor());
00455 #else
00456     PDBAtomPtr *atomptr = atomArray;
00457     PDBAtom *atom = *atomptr;
00458     Vector p(atom->xcoor(),atom->ycoor(),atom->zcoor());
00459 #endif
00460     Vector s(lattice.a_r()*p, lattice.b_r()*p, lattice.c_r()*p);
00461     smin = smax = s;
00462     atomptr++;
00463     for(int i=1; i<atomCount; i++, atomptr++){
00464 #ifdef MEM_OPT_VERSION
00465       p = Vector(atomptr->xcoor(),atomptr->ycoor(),atomptr->zcoor());
00466 #else
00467       atom = *atomptr;
00468       p = Vector (atom->xcoor(),atom->ycoor(),atom->zcoor());
00469 #endif
00470       s = Vector(lattice.a_r()*p, lattice.b_r()*p, lattice.c_r()*p);
00471       if      (smin.x > s.x) smin.x = s.x;
00472       else if (smax.x < s.x) smax.x = s.x;
00473       if      (smin.y > s.y) smin.y = s.y;
00474       else if (smax.y < s.y) smax.y = s.y;
00475       if      (smin.z > s.z) smin.z = s.z;
00476       else if (smax.z < s.z) smax.z = s.z;
00477     }
00478   }
00479   // shift using the origin
00480   BigReal origin_shift;
00481   origin_shift = lattice.a_r() * lattice.origin();
00482   smin.x -= origin_shift;
00483   smax.x -= origin_shift;
00484   origin_shift = lattice.b_r() * lattice.origin();
00485   smin.y -= origin_shift;
00486   smax.y -= origin_shift;
00487   origin_shift = lattice.c_r() * lattice.origin();
00488   smin.z -= origin_shift;
00489   smax.z -= origin_shift;
00490 }
00491 
00492 //#define TEST_PDB_CLASS
00493 #ifdef TEST_PDB_CLASS
00494 
00495 main()
00496 {
00497  PDB *pdb = new PDB("pti.pdb");
00498  if ( atomArray == NULL )
00499  {
00500    NAMD_die("memory allocation failed in main of test PDB class");
00501  }
00502  ilist = pdb->find_atom_name("CA");
00503  if (!ilist)
00504    printf("None found.\n");
00505  else {
00506    int i;
00507    char s[200];
00508    for (i=0; i<ilist->num(); i++) {
00509      pdb->atom((*ilist)[i]) -> sprint(s);
00510      printf("%s\n", s);
00511    }
00512    delete ilist;
00513  } // if ilist
00514  
00515  printf("Now, search through space.\n");
00516  
00517  ilist = pdb->find_atoms_in_region(4.38, 19.5, 3.0,
00518                                    4.40, 20.0, 3.2 );
00519  if (!ilist)
00520    printf("None found.\n");
00521  else {
00522    int i;
00523    char s[200];
00524    printf("%d found\n", ilist -> num());
00525    for (i=0; i<ilist->num(); i++) {
00526      pdb->atom((*ilist)[i]) -> sprint(s);
00527      printf("%s\n", s);
00528    }
00529    delete ilist;
00530  }
00531 }
00532 #endif // TEST_PDB_CLASS
00533 
00534 
00535 
00536 // This function was borrowed from VMD code in "ReadPARM.C".
00537 // It skips to a new line.
00538 static int readtoeoln(FILE *f) {
00539   char c;
00540 
00541   /* skip to eoln */
00542   while((c = getc(f)) != '\n') {
00543     if (c == EOF) 
00544       return -1;
00545   }
00546 
00547   return 0;
00548 }  
00549 
00550 
00551 
00552 // read in an AMBER coordinate file and populate the PDB structure
00553 PDB::PDB( const char *filename, Ambertoppar *amber_data)
00554 { int i,j,k;
00555   BigReal coor[3];
00556   char buf[13],resname[5],atomname[5];
00557   FILE *infile;
00558   PDBAtom *pdb;
00559 
00560   if ((infile=Fopen(filename, "r")) == NULL)
00561     NAMD_err("Can't open AMBER coordinate file!");
00562 
00563   readtoeoln(infile);  // Skip the first line (title)
00564 
00565   fscanf(infile,"%d",&atomCount);  // Read in num of atoms
00566   if (atomCount != amber_data->Natom)
00567     NAMD_die("Num of atoms in coordinate file is different from that in parm file!");
00568   readtoeoln(infile);
00569 
00570 #ifdef MEM_OPT_VERSION
00571   altlocArray = 0;
00572   atomArray = new PDBCoreData[atomCount];
00573 #else
00574   atomArray = new PDBAtomPtr[atomCount];
00575 #endif
00576 
00577   if ( atomArray == NULL )
00578   {
00579     NAMD_die("memory allocation failed in PDB::PDB");
00580   }
00581   
00582   // Read in the coordinates, which are in the format of 6F12.7
00583   // Other fields are copied from "amber_data"
00584   for (i=0; i<atomCount; ++i)
00585   { // Read x,y,z coordinates
00586     for (j=0; j<3; ++j)
00587     { for (k=0; k<12; ++k)
00588       { buf[k]=getc(infile);
00589         if (buf[k]=='\n' || buf[k]=='\0' || buf[k]==EOF)
00590           NAMD_die("Error reading AMBER coordinate file!");
00591       }
00592       buf[12] = '\0';
00593       coor[j] = atof(buf);
00594     }
00595     if (i%2 == 1)
00596       readtoeoln(infile);
00597 #ifdef MEM_OPT_VERSION
00598     atomArray[i].coor[0] = coor[0];
00599     atomArray[i].coor[1] = coor[1];
00600     atomArray[i].coor[2] = coor[2];
00601     atomArray[i].myoccupancy = PDBAtom::default_occupancy;
00602     atomArray[i].tempfactor = PDBAtom::default_temperaturefactor;
00603 #else
00604     // Copy name, resname and resid from "amber_data"
00605     for (j=0; j<4; ++j)
00606     { resname[j] = amber_data->ResNames[amber_data->AtomRes[i]*4+j];
00607       atomname[j] = amber_data->AtomNames[i*4+j];
00608     }
00609     resname[4] = atomname[4] = '\0';
00610     // Create a new PDB record, and fill in its entries
00611     pdb = new PDBAtomRecord("");
00612     pdb->name(atomname);
00613     pdb->residuename(resname);
00614     pdb->serialnumber(i+1);
00615     pdb->residueseq(amber_data->AtomRes[i]+1);
00616     pdb->coordinates(coor);
00617     atomArray[i] = pdb;  // Include the new record into the array
00618 #endif
00619   }
00620 }
00621 
00622 #define LINESIZE 100
00623 
00624 /* This constructor initializes the PDB data using a Gromacs
00625    coordinate file, generating an error message if the file
00626    can't be parsed or if its contents don't jive with what is in
00627    the topo file <topology>. */
00628 PDB::PDB(const char *filename, const GromacsTopFile *topology) {
00629   int i;
00630   char buf[LINESIZE];
00631   FILE *infile;
00632   
00633   /* open up the coordinate file */
00634   infile=Fopen(filename, "r");
00635   if (infile == NULL)
00636     NAMD_err("Can't open GROMACS coordinate file!");
00637 
00638   fgets(buf,LINESIZE-1,infile); // get the title
00639   // if(strcmp(buf,topology->getSystemName()) != 0)
00640   //   NAMD_die("System names in topology and coordinate files differ.");
00641 
00642   fgets(buf,LINESIZE-1,infile); // get the number of atoms
00643   sscanf(buf,"%d",&atomCount);
00644   if (atomCount != topology->getNumAtoms())
00645     NAMD_die("Num of atoms in coordinate file is different from that in topology file!");
00646 
00647   /* read in the atoms */
00648 #ifdef MEM_OPT_VERSION
00649   altlocArray = 0;
00650   atomArray = new PDBCoreData[atomCount];
00651 #else
00652   atomArray = new PDBAtomPtr[atomCount];
00653 #endif
00654   if ( atomArray == NULL )
00655     NAMD_die("memory allocation failed in PDB::PDB");
00656 
00657 #ifdef MEM_OPT_VERSION
00658   for(i=0; i<atomCount; i++){
00659       fgets(buf,LINESIZE-1,infile); // get a line
00660       char *buf2 = buf+20; // skip three fields to get to the coordinates
00661       BigReal coor[3];
00662       if(3 != sscanf(buf2,"%lf%lf%lf", &coor[0],&coor[1],&coor[2]))
00663           NAMD_die("Couldn't get three coordinates from file.");
00664 
00665       coor[0] *= 10; // convert to angstroms from nanometers
00666       coor[1] *= 10;
00667       coor[2] *= 10;
00668 
00669       atomArray[i].coor[0] = coor[0];
00670       atomArray[i].coor[1] = coor[1];
00671       atomArray[i].coor[2] = coor[2];
00672       atomArray[i].myoccupancy = PDBAtom::default_occupancy;
00673       atomArray[i].tempfactor = PDBAtom::default_temperaturefactor;
00674   }
00675 #else
00676   for (i=0;i<atomCount;i++) {
00677     char *buf2, resname[11], atomname[11], atmtype[11];
00678     int resnum, typenum;
00679     Real charge,mass;
00680     BigReal coor[3];
00681     PDBAtom *pdb = new PDBAtomRecord("");  
00682     
00683     fgets(buf,LINESIZE-1,infile); // get a line
00684     buf2 = buf+20; // skip three fields to get to the coordinates
00685     if(3 != sscanf(buf2,"%lf%lf%lf",
00686                    &coor[0],&coor[1],&coor[2]))
00687       NAMD_die("Couldn't get three coordinates from file.");
00688     topology->getAtom(i,&resnum,resname,
00689                       atomname,atmtype,&typenum,&charge,&mass);
00690     coor[0] *= 10; // convert to angstroms from nanometers
00691     coor[1] *= 10;
00692     coor[2] *= 10;
00693     
00694     pdb->name(atomname);
00695     pdb->residuename(resname);
00696     pdb->serialnumber(i+1);
00697     pdb->residueseq(resnum+1);
00698     pdb->coordinates(coor);
00699     
00700     atomArray[i] = pdb;  // Include the new record into the array
00701   }
00702 #endif
00703 }
00704 
00705 #ifdef MEM_OPT_VERSION
00706 void PDBCoreData::sprint( char *outstr, PDBData::PDBFormatStyle usestyle){
00707     if(usestyle == PDBData::COLUMNS){
00708         //sprint_columns copied from PDBAtom::sprint_columns
00709         for(int i=0; i<79; i++) outstr[i] = 32;
00710         outstr[79] = 0;
00711         PDBData::sprintcol(outstr, PDBAtom::SX, PDBAtom::LCOOR, PDBAtom::LCOORPREC, xcoor());
00712         PDBData::sprintcol(outstr, PDBAtom::SY, PDBAtom::LCOOR, PDBAtom::LCOORPREC, ycoor());
00713         PDBData::sprintcol(outstr, PDBAtom::SZ, PDBAtom::LCOOR, PDBAtom::LCOORPREC, zcoor());
00714         PDBData::sprintcol(outstr, PDBAtom::SOCC, PDBAtom::LOCC, PDBAtom::LOCCPREC, occupancy());
00715         PDBData::sprintcol(outstr, PDBAtom::STEMPF, PDBAtom::LTEMPF, PDBAtom::LTEMPFPREC, temperaturefactor());
00716     }else{
00717         //sprint_fields
00718         char tmpstr[50];        
00719         if (xcoor() == PDBAtom::default_coor)
00720             sprintf(tmpstr, " #");
00721            else
00722             sprintf(tmpstr, " %*.*f", PDBAtom::LCOOR, PDBAtom::LCOORPREC, xcoor());
00723         strcat(outstr, tmpstr);
00724         if (ycoor() == PDBAtom::default_coor)
00725             sprintf(tmpstr, " #");
00726            else
00727             sprintf(tmpstr, " %*.*f", PDBAtom::LCOOR, PDBAtom::LCOORPREC, ycoor());
00728         strcat(outstr, tmpstr);
00729         if (zcoor() == PDBAtom::default_coor)
00730             sprintf(tmpstr, " #");
00731            else
00732             sprintf(tmpstr, " %*.*f", PDBAtom::LCOOR, PDBAtom::LCOORPREC, zcoor());
00733         strcat(outstr, tmpstr);
00734       
00735         sprintf(tmpstr, " %*.*f",  PDBAtom::LOCC, PDBAtom::LOCCPREC, occupancy());
00736         strcat(outstr, tmpstr);
00737       
00738         sprintf(tmpstr, " %*.*f", PDBAtom::LTEMPF, PDBAtom::LTEMPFPREC, temperaturefactor());
00739         strcat(outstr, tmpstr);        
00740     }
00741 }
00742 #endif
00743 

Generated on Mon Nov 20 01:17:14 2017 for NAMD by  doxygen 1.4.7