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 = atomAlloc = 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   atomAlloc = 0;
00218   atomArray = new PDBAtomPtr[atomCount];
00219   if ( atomArray == NULL )
00220   {
00221     NAMD_die("memory allocation failed in PDB::PDB");
00222   }
00223   PDBAtomList *tmp = atomListHead;
00224   int i=0;                              // just need to copy the pointers
00225   for (i=0, tmp = atomListHead; tmp != NULL; tmp = tmp -> next, i++) {
00226     atomArray[i] = tmp -> data;
00227   }
00228   // now delete the linked list (w/o deleting the data)
00229   PDBAtomList *tmp2;
00230   for (tmp2 = tmp = atomListHead; tmp != NULL; tmp = tmp2) {
00231     tmp2 = tmp->next;
00232     delete tmp;
00233   }
00234   atomListHead = atomListTail = NULL;
00235 #endif 
00236  }  // everything converted
00237  
00238 }
00239 
00240 //  Destructor - delete all the data pointed to by the array
00241 //   and then delete the array
00242 PDB::~PDB( void )
00243 {
00244 #ifndef MEM_OPT_VERSION
00245         int i;
00246         if ( atomAlloc ) {
00247           delete [] atomAlloc;
00248         } else {
00249           for (i=atomCount-1; i>=0; i--)
00250             delete atomArray[i];
00251         }
00252 #else
00253         delete [] altlocArray;
00254 #endif
00255         delete [] atomArray;
00256         atomArray = NULL;
00257         atomCount = 0;
00258 }
00259 
00260 // print the PDB file out to a given file name
00261 void PDB::write(const char *outfilename, const char *commentline)
00262 {
00263         int i;
00264         char s[200];
00265         FILE *outfile;
00266         if ((outfile = fopen(outfilename, "w")) == NULL) {
00267            sprintf(s, "Cannot open file '%s' in PDB::write.", outfilename);
00268            NAMD_err(s);
00269         }
00270 
00271         if (commentline != NULL)
00272         {
00273                 sprintf(s, "REMARK  %s\n", commentline);
00274                 if (fputs(s, outfile) == EOF)
00275                 {
00276                         NAMD_err("EOF in PDB::write writing the comment line - file system full?");
00277                 }
00278         }
00279 
00280         for (i=0; i<atomCount; i++){ // I only contain ATOM/HETATM records
00281 #ifdef MEM_OPT_VERSION    
00282       atomArray[i].sprint(s, PDBData::COLUMNS);  
00283 #else
00284           atomArray[i]->sprint(s, PDBData::COLUMNS);
00285 #endif
00286           if ( (fputs(s, outfile)    == EOF) || 
00287                (fputs("\n", outfile) == EOF)    ) {
00288             sprintf(s, "EOF in PDB::write line %d - file system full?", i);
00289             NAMD_err(s);
00290           }
00291         }
00292         if (fputs("END\n", outfile) == EOF) {
00293            NAMD_err("EOF in PDB::write while printing 'END' -- file system full?");
00294         }
00295         if (fclose(outfile) == EOF) {
00296            NAMD_err("EOF in PDB::write while closing -- file system full?");
00297         }
00298           
00299 }
00300 
00301 // store the info on the linked list
00302 void PDB::add_atom_element( PDBAtom *newAtom)
00303 {
00304   PDBAtomList *tmp = new PDBAtomList;
00305   if ( tmp == NULL )
00306   {
00307     NAMD_die("memory allocation failed in PDB::add_atom_element");
00308   }
00309   tmp -> data = newAtom;
00310   tmp -> next = NULL;
00311   
00312   if (atomListHead == NULL) {        // make the list
00313     atomListHead = atomListTail = tmp;
00314   } else {
00315     atomListTail -> next = tmp;       // add to the tail
00316     atomListTail = tmp;
00317   }
00318   atomCount++;
00319 }
00320 
00321 
00322 // return the number of atoms found
00323 int PDB::num_atoms( void)
00324 {
00325   return atomCount;
00326 }
00327 
00328 
00329 // Reset all the atom positions.  This is used in preparation for
00330 // output in cases like the restart files, etc.
00331 void PDB::set_all_positions(Vector *pos)
00332 {
00333 #ifdef MEM_OPT_VERSION
00334     for(int i=0; i<atomCount; i++){
00335         atomArray[i].coor[0] = pos[i].x;
00336         atomArray[i].coor[1] = pos[i].y;
00337         atomArray[i].coor[2] = pos[i].z;
00338     }
00339 #else
00340         int i;
00341         PDBAtomPtr *atomptr;
00342 
00343         for (i=0, atomptr=atomArray; i<atomCount; atomptr++, i++)
00344         {
00345                 (*atomptr)->xcoor(pos[i].x);
00346                 (*atomptr)->ycoor(pos[i].y);
00347                 (*atomptr)->zcoor(pos[i].z);
00348         }
00349 #endif
00350 }
00351 
00352 void PDB::get_position_for_atom(Vector *pos, int aid){
00353 #ifdef MEM_OPT_VERSION
00354     pos->x = atomArray[aid].coor[0];
00355     pos->y = atomArray[aid].coor[1];
00356     pos->z = atomArray[aid].coor[2];
00357 #else
00358     PDBAtomPtr *atomptr = &atomArray[aid];
00359     pos->x = (*atomptr)->xcoor();
00360     pos->y = (*atomptr)->ycoor();
00361     pos->z = (*atomptr)->zcoor();
00362 #endif
00363 }
00364 
00365 //  Get all the atom positions into a list of Vectors
00366 void PDB::get_all_positions(Vector *pos)
00367 {
00368 #ifdef MEM_OPT_VERSION
00369     for(int i=0; i<atomCount; i++){
00370         pos[i].x = atomArray[i].coor[0];
00371         pos[i].y = atomArray[i].coor[1];
00372         pos[i].z = atomArray[i].coor[2];        
00373     }
00374 #else
00375         int i;
00376         PDBAtomPtr *atomptr;
00377 
00378         for (i=0, atomptr=atomArray; i<atomCount; atomptr++, i++)
00379         {
00380                 pos[i].x = (*atomptr)->xcoor();
00381                 pos[i].y = (*atomptr)->ycoor();
00382                 pos[i].z = (*atomptr)->zcoor();
00383         }
00384 #endif
00385 }
00386 
00387 //  given an index, return that atom
00388 #ifdef MEM_OPT_VERSION
00389 PDBCoreData *PDB::atom(int place){
00390     if(place<0 || place>=atomCount) return NULL;
00391     return &atomArray[place];
00392 }
00393 #else
00394 PDBAtom *PDB::atom(int place)
00395 {
00396   if (place <0 || place >= atomCount)
00397     return NULL;
00398   return atomArray[place];
00399 }
00400 #endif
00401 
00402 
00403 // find the lowest and highest bounds based on a fraction of the atoms
00404 //void PDB::find_extremes(BigReal *min, BigReal *max, Vector rec, BigReal frac) const
00405 void PDB::find_extremes_helper(
00406     SortableResizeArray<BigReal> &coor,
00407     BigReal &min, BigReal &max, Vector rec, BigReal frac
00408     )
00409 {
00410     SortableResizeArray<BigReal>::iterator c_i = coor.begin();
00411 #ifdef MEM_OPT_VERSION
00412     PDBCoreData *atomptr = atomArray;
00413     for(int i=0; i<atomCount; i++, atomptr++){
00414         Vector pos(atomptr->xcoor(),atomptr->ycoor(),atomptr->zcoor());
00415         c_i[i] = rec*pos;
00416     }
00417 #else
00418     PDBAtomPtr *atomptr = atomArray;
00419     for (int i=0; i<atomCount; ++i, ++atomptr) {
00420       PDBAtom *atom = *atomptr;
00421       Vector pos(atom->xcoor(),atom->ycoor(),atom->zcoor());
00422       c_i[i] = rec*pos;
00423     }
00424 #endif
00425     coor.sort();
00426     int ilow = (int)((1.0 - frac) * atomCount);
00427     if ( ilow < 0 ) ilow = 0;
00428     if ( ilow > atomCount/2 ) ilow = atomCount/2;
00429     int ihigh = atomCount - ilow - 1;
00430     BigReal span = coor[ihigh] - coor[ilow];
00431     BigReal extension = (1.0 - frac) * span / (2.0 * frac - 1.0);
00432     max = coor[ihigh] + extension;
00433     min = coor[ilow] - extension;
00434 }
00435 
00436 // Find the extreme edges of molecule in scaled coordinates,
00437 // where "frac" sets bounds based on a fraction of the atoms.
00438 void PDB::find_extremes(const Lattice &lattice, BigReal frac) {
00439   if (atomCount == 0) {
00440     smin = smax = 0;
00441   }
00442   else if (frac < 1.0) {
00443     // for now use the previous "sort the array" approach
00444     // for solving "frac"th largest and smallest selection problems
00445     SortableResizeArray<BigReal> coor;
00446     coor.resize(atomCount);  // allocate array space once
00447     find_extremes_helper(coor, smin.x, smax.x, lattice.a_r(), frac);
00448     find_extremes_helper(coor, smin.y, smax.y, lattice.b_r(), frac);
00449     find_extremes_helper(coor, smin.z, smax.z, lattice.c_r(), frac);
00450   }
00451   else {
00452     // finding absolute min and max does not require sorting
00453 #ifdef MEM_OPT_VERSION
00454     PDBCoreData *atomptr = atomArray;
00455     Vector p(atomptr->xcoor(),atomptr->ycoor(),atomptr->zcoor());
00456 #else
00457     PDBAtomPtr *atomptr = atomArray;
00458     PDBAtom *atom = *atomptr;
00459     Vector p(atom->xcoor(),atom->ycoor(),atom->zcoor());
00460 #endif
00461     Vector s(lattice.a_r()*p, lattice.b_r()*p, lattice.c_r()*p);
00462     smin = smax = s;
00463     atomptr++;
00464     for(int i=1; i<atomCount; i++, atomptr++){
00465 #ifdef MEM_OPT_VERSION
00466       p = Vector(atomptr->xcoor(),atomptr->ycoor(),atomptr->zcoor());
00467 #else
00468       atom = *atomptr;
00469       p = Vector (atom->xcoor(),atom->ycoor(),atom->zcoor());
00470 #endif
00471       s = Vector(lattice.a_r()*p, lattice.b_r()*p, lattice.c_r()*p);
00472       if      (smin.x > s.x) smin.x = s.x;
00473       else if (smax.x < s.x) smax.x = s.x;
00474       if      (smin.y > s.y) smin.y = s.y;
00475       else if (smax.y < s.y) smax.y = s.y;
00476       if      (smin.z > s.z) smin.z = s.z;
00477       else if (smax.z < s.z) smax.z = s.z;
00478     }
00479   }
00480   // shift using the origin
00481   BigReal origin_shift;
00482   origin_shift = lattice.a_r() * lattice.origin();
00483   smin.x -= origin_shift;
00484   smax.x -= origin_shift;
00485   origin_shift = lattice.b_r() * lattice.origin();
00486   smin.y -= origin_shift;
00487   smax.y -= origin_shift;
00488   origin_shift = lattice.c_r() * lattice.origin();
00489   smin.z -= origin_shift;
00490   smax.z -= origin_shift;
00491 }
00492 
00493 //#define TEST_PDB_CLASS
00494 #ifdef TEST_PDB_CLASS
00495 
00496 main()
00497 {
00498  PDB *pdb = new PDB("pti.pdb");
00499  if ( atomArray == NULL )
00500  {
00501    NAMD_die("memory allocation failed in main of test PDB class");
00502  }
00503  ilist = pdb->find_atom_name("CA");
00504  if (!ilist)
00505    printf("None found.\n");
00506  else {
00507    int i;
00508    char s[200];
00509    for (i=0; i<ilist->num(); i++) {
00510      pdb->atom((*ilist)[i]) -> sprint(s);
00511      printf("%s\n", s);
00512    }
00513    delete ilist;
00514  } // if ilist
00515  
00516  printf("Now, search through space.\n");
00517  
00518  ilist = pdb->find_atoms_in_region(4.38, 19.5, 3.0,
00519                                    4.40, 20.0, 3.2 );
00520  if (!ilist)
00521    printf("None found.\n");
00522  else {
00523    int i;
00524    char s[200];
00525    printf("%d found\n", ilist -> num());
00526    for (i=0; i<ilist->num(); i++) {
00527      pdb->atom((*ilist)[i]) -> sprint(s);
00528      printf("%s\n", s);
00529    }
00530    delete ilist;
00531  }
00532 }
00533 #endif // TEST_PDB_CLASS
00534 
00535 
00536 
00537 // This function was borrowed from VMD code in "ReadPARM.C".
00538 // It skips to a new line.
00539 static int readtoeoln(FILE *f) {
00540   char c;
00541 
00542   /* skip to eoln */
00543   while((c = getc(f)) != '\n') {
00544     if (c == EOF) 
00545       return -1;
00546   }
00547 
00548   return 0;
00549 }  
00550 
00551 
00552 
00553 // read in an AMBER coordinate file and populate the PDB structure
00554 PDB::PDB( const char *filename, Ambertoppar *amber_data)
00555 { int i,j,k;
00556   BigReal coor[3];
00557   char buf[13],resname[5],atomname[5];
00558   FILE *infile;
00559   PDBAtom *pdb;
00560 
00561   if ((infile=Fopen(filename, "r")) == NULL)
00562     NAMD_err("Can't open AMBER coordinate file!");
00563 
00564   readtoeoln(infile);  // Skip the first line (title)
00565 
00566   fscanf(infile,"%d",&atomCount);  // Read in num of atoms
00567   if (atomCount != amber_data->Natom)
00568     NAMD_die("Num of atoms in coordinate file is different from that in parm file!");
00569   readtoeoln(infile);
00570 
00571 #ifdef MEM_OPT_VERSION
00572   altlocArray = 0;
00573   atomArray = new PDBCoreData[atomCount];
00574 #else
00575   atomAlloc = 0;
00576   atomArray = new PDBAtomPtr[atomCount];
00577 #endif
00578 
00579   if ( atomArray == NULL )
00580   {
00581     NAMD_die("memory allocation failed in PDB::PDB");
00582   }
00583   
00584   // Read in the coordinates, which are in the format of 6F12.7
00585   // Other fields are copied from "amber_data"
00586   for (i=0; i<atomCount; ++i)
00587   { // Read x,y,z coordinates
00588     for (j=0; j<3; ++j)
00589     { for (k=0; k<12; ++k)
00590       { buf[k]=getc(infile);
00591         if (buf[k]=='\n' || buf[k]=='\0' || buf[k]==EOF)
00592           NAMD_die("Error reading AMBER coordinate file!");
00593       }
00594       buf[12] = '\0';
00595       coor[j] = atof(buf);
00596     }
00597     if (i%2 == 1)
00598       readtoeoln(infile);
00599 #ifdef MEM_OPT_VERSION
00600     atomArray[i].coor[0] = coor[0];
00601     atomArray[i].coor[1] = coor[1];
00602     atomArray[i].coor[2] = coor[2];
00603     atomArray[i].myoccupancy = PDBAtom::default_occupancy;
00604     atomArray[i].tempfactor = PDBAtom::default_temperaturefactor;
00605 #else
00606     // Copy name, resname and resid from "amber_data"
00607     for (j=0; j<4; ++j)
00608     { resname[j] = amber_data->ResNames[amber_data->AtomRes[i]*4+j];
00609       atomname[j] = amber_data->AtomNames[i*4+j];
00610     }
00611     resname[4] = atomname[4] = '\0';
00612     // Create a new PDB record, and fill in its entries
00613     pdb = new PDBAtomRecord("");
00614     pdb->name(atomname);
00615     pdb->residuename(resname);
00616     pdb->serialnumber(i+1);
00617     pdb->residueseq(amber_data->AtomRes[i]+1);
00618     pdb->coordinates(coor);
00619     atomArray[i] = pdb;  // Include the new record into the array
00620 #endif
00621   }
00622 }
00623 
00624 #define LINESIZE 100
00625 
00626 /* This constructor initializes the PDB data using a Gromacs
00627    coordinate file, generating an error message if the file
00628    can't be parsed or if its contents don't jive with what is in
00629    the topo file <topology>. */
00630 PDB::PDB(const char *filename, const GromacsTopFile *topology) {
00631   int i;
00632   char buf[LINESIZE];
00633   FILE *infile;
00634   
00635   /* open up the coordinate file */
00636   infile=Fopen(filename, "r");
00637   if (infile == NULL)
00638     NAMD_err("Can't open GROMACS coordinate file!");
00639 
00640   fgets(buf,LINESIZE-1,infile); // get the title
00641   // if(strcmp(buf,topology->getSystemName()) != 0)
00642   //   NAMD_die("System names in topology and coordinate files differ.");
00643 
00644   fgets(buf,LINESIZE-1,infile); // get the number of atoms
00645   sscanf(buf,"%d",&atomCount);
00646   if (atomCount != topology->getNumAtoms())
00647     NAMD_die("Num of atoms in coordinate file is different from that in topology file!");
00648 
00649   /* read in the atoms */
00650 #ifdef MEM_OPT_VERSION
00651   altlocArray = 0;
00652   atomArray = new PDBCoreData[atomCount];
00653 #else
00654   atomAlloc = 0;
00655   atomArray = new PDBAtomPtr[atomCount];
00656 #endif
00657   if ( atomArray == NULL )
00658     NAMD_die("memory allocation failed in PDB::PDB");
00659 
00660 #ifdef MEM_OPT_VERSION
00661   for(i=0; i<atomCount; i++){
00662       fgets(buf,LINESIZE-1,infile); // get a line
00663       char *buf2 = buf+20; // skip three fields to get to the coordinates
00664       BigReal coor[3];
00665       if(3 != sscanf(buf2,"%lf%lf%lf", &coor[0],&coor[1],&coor[2]))
00666           NAMD_die("Couldn't get three coordinates from file.");
00667 
00668       coor[0] *= 10; // convert to angstroms from nanometers
00669       coor[1] *= 10;
00670       coor[2] *= 10;
00671 
00672       atomArray[i].coor[0] = coor[0];
00673       atomArray[i].coor[1] = coor[1];
00674       atomArray[i].coor[2] = coor[2];
00675       atomArray[i].myoccupancy = PDBAtom::default_occupancy;
00676       atomArray[i].tempfactor = PDBAtom::default_temperaturefactor;
00677   }
00678 #else
00679   for (i=0;i<atomCount;i++) {
00680     char *buf2, resname[11], atomname[11], atmtype[11];
00681     int resnum, typenum;
00682     Real charge,mass;
00683     BigReal coor[3];
00684     PDBAtom *pdb = new PDBAtomRecord("");  
00685     
00686     fgets(buf,LINESIZE-1,infile); // get a line
00687     buf2 = buf+20; // skip three fields to get to the coordinates
00688     if(3 != sscanf(buf2,"%lf%lf%lf",
00689                    &coor[0],&coor[1],&coor[2]))
00690       NAMD_die("Couldn't get three coordinates from file.");
00691     topology->getAtom(i,&resnum,resname,
00692                       atomname,atmtype,&typenum,&charge,&mass);
00693     coor[0] *= 10; // convert to angstroms from nanometers
00694     coor[1] *= 10;
00695     coor[2] *= 10;
00696     
00697     pdb->name(atomname);
00698     pdb->residuename(resname);
00699     pdb->serialnumber(i+1);
00700     pdb->residueseq(resnum+1);
00701     pdb->coordinates(coor);
00702     
00703     atomArray[i] = pdb;  // Include the new record into the array
00704   }
00705 #endif
00706 }
00707 
00708 #ifdef MEM_OPT_VERSION
00709 void PDBCoreData::sprint( char *outstr, PDBData::PDBFormatStyle usestyle){
00710     if(usestyle == PDBData::COLUMNS){
00711         //sprint_columns copied from PDBAtom::sprint_columns
00712         for(int i=0; i<79; i++) outstr[i] = 32;
00713         outstr[79] = 0;
00714         PDBData::sprintcol(outstr, PDBAtom::SX, PDBAtom::LCOOR, PDBAtom::LCOORPREC, xcoor());
00715         PDBData::sprintcol(outstr, PDBAtom::SY, PDBAtom::LCOOR, PDBAtom::LCOORPREC, ycoor());
00716         PDBData::sprintcol(outstr, PDBAtom::SZ, PDBAtom::LCOOR, PDBAtom::LCOORPREC, zcoor());
00717         PDBData::sprintcol(outstr, PDBAtom::SOCC, PDBAtom::LOCC, PDBAtom::LOCCPREC, occupancy());
00718         PDBData::sprintcol(outstr, PDBAtom::STEMPF, PDBAtom::LTEMPF, PDBAtom::LTEMPFPREC, temperaturefactor());
00719     }else{
00720         //sprint_fields
00721         char tmpstr[50];        
00722         if (xcoor() == PDBAtom::default_coor)
00723             sprintf(tmpstr, " #");
00724            else
00725             sprintf(tmpstr, " %*.*f", PDBAtom::LCOOR, PDBAtom::LCOORPREC, xcoor());
00726         strcat(outstr, tmpstr);
00727         if (ycoor() == PDBAtom::default_coor)
00728             sprintf(tmpstr, " #");
00729            else
00730             sprintf(tmpstr, " %*.*f", PDBAtom::LCOOR, PDBAtom::LCOORPREC, ycoor());
00731         strcat(outstr, tmpstr);
00732         if (zcoor() == PDBAtom::default_coor)
00733             sprintf(tmpstr, " #");
00734            else
00735             sprintf(tmpstr, " %*.*f", PDBAtom::LCOOR, PDBAtom::LCOORPREC, zcoor());
00736         strcat(outstr, tmpstr);
00737       
00738         sprintf(tmpstr, " %*.*f",  PDBAtom::LOCC, PDBAtom::LOCCPREC, occupancy());
00739         strcat(outstr, tmpstr);
00740       
00741         sprintf(tmpstr, " %*.*f", PDBAtom::LTEMPF, PDBAtom::LTEMPFPREC, temperaturefactor());
00742         strcat(outstr, tmpstr);        
00743     }
00744 }
00745 #endif
00746 

Generated on Fri Jul 20 01:17:15 2018 for NAMD by  doxygen 1.4.7