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

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

Generated on Sat May 25 04:07:19 2013 for NAMD by  doxygen 1.3.9.1