EMAN2
pdbreader.cpp
Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Muthu Alagappan, 08/11/2004 (m.alagappan901@gmail.com)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #include "pdbreader.h"
00037 #include <vector>
00038 #include <cstring>
00039 #include <string>
00040 #include <stdio.h>
00041 #include <iostream>
00042 
00043 using namespace EMAN;
00044 
00045 
00046 PDBReader::PDBReader()
00047 {
00048         points = 0;
00049         n = 0;
00050 }
00051 
00052 PDBReader::PDBReader( int nn)
00053 {
00054         n = nn;
00055         points = (double *) calloc(4 * n, sizeof(double));
00056         ter_stop = 0;
00057 }
00058 
00059 PDBReader::~PDBReader()
00060 {
00061         if( points )
00062         {
00063                 free(points);
00064                 points = 0;
00065         }
00066 }
00067 
00068 void PDBReader::zero()
00069 {
00070         memset((void *) points, 0, 4 * n * sizeof(double));
00071 }
00072 
00073 
00074 //copies all relevant information that exists into a new PDBReader object.
00075 PDBReader *PDBReader::copy()
00076 {
00077         PDBReader *pa2 = new PDBReader();
00078         pa2->set_number_points(get_number_points());
00079         double *pa2data = pa2->get_points_array();
00080         memcpy(pa2data, get_points_array(), sizeof(double) * 4 * get_number_points());
00081         pa2->pWords = pWords;
00082         pa2->atomName = atomName;
00083         pa2->residueName = residueName;
00084         pa2->chainId = chainId;
00085         pa2->elementSym = elementSym;
00086         pa2->tail = tail;
00087         pa2->head = head;
00088         pa2->pointInfo = pointInfo;
00089         pa2->lines = lines;
00090         return pa2;
00091 }
00092 
00093 
00094 PDBReader & PDBReader::operator=(PDBReader & pa)
00095 {
00096         if (this != &pa) {
00097                 set_number_points(pa.get_number_points());
00098                 memcpy(get_points_array(), pa.get_points_array(), sizeof(double) * 4 * get_number_points());
00099         }
00100         return *this;
00101 }
00102 
00103 size_t PDBReader::get_number_points() const
00104 {
00105         return n;
00106 }
00107 
00108 //reallocates the number of points
00109 void PDBReader::set_number_points(size_t nn)
00110 {
00111         if (n != nn) {
00112                 n = nn;
00113                 points = (double *) realloc(points, 4 * n * sizeof(double));
00114         }
00115 }
00116 
00117 double *PDBReader::get_points_array()
00118 {
00119         return points;
00120 }
00121 
00122 void PDBReader::set_points_array(double *p)
00123 {
00124         points = p;
00125 }
00126 
00127 
00128 
00129 // adds purely the x values to a vector
00130 vector<float> PDBReader::get_x() {
00131         if (count_stop == 0) {count_stop = atomName.size();}
00132         for (int i=0; i<count_stop; i++) {
00133                 x.push_back((float)points[4*i]);
00134                 y.push_back((float)points[4*i + 1]);
00135                 z.push_back((float)points[4*i + 2]);
00136                 resNum.push_back(pointInfo[2*i+1]);
00137         }
00138 
00139         return x;
00140 }
00141 
00142 // adds purely the y values to a vector
00143 vector<float> PDBReader::get_y() {
00144         //cout << y.size() << endl;
00145         return y;
00146 }
00147 
00148 // adds purely the z values to a vector
00149 vector<float> PDBReader::get_z() {
00150         return z;
00151 }
00152 
00153 vector<string> PDBReader::get_atomName() {
00154         return atomName;
00155 }
00156 
00157 vector<string> PDBReader::get_resName() {
00158         return residueName;
00159 }
00160 
00161 vector<int> PDBReader::get_resNum() {
00162         return resNum;
00163 }
00164 
00165 
00166 //Accurately parses a pdb file for all information
00167 bool PDBReader::read_from_pdb(const char *file)
00168 {
00169 
00170         pWords.clear();
00171         atomName.clear();
00172         residueName.clear();
00173         chainId.clear();
00174         elementSym.clear();
00175         tail.clear();
00176         head.clear();
00177         lines.clear();
00178         pointInfo.clear();
00179         ter_stop = 0;
00180         count_stop = 0;
00181 
00182         struct stat filestat;
00183         stat(file, &filestat);
00184         set_number_points(( int)(filestat.st_size / 80 + 1)); //80 bytes per line
00185 #ifdef DEBUG
00186         printf("PDBReader::read_pdb(): try %4lu atoms first\n", get_number_points());
00187 #endif
00188 
00189         FILE *fp = fopen(file, "r");
00190         if(!fp) {
00191                 fprintf(stderr,"ERROR in PDBReader::read_pdb(): cannot open file %s\n",file);
00192                 throw;
00193         }
00194         char s[200];
00195         size_t count = 0;
00196 
00197         while ((fgets(s, 200, fp) != NULL)) {
00198                 lines.push_back(s);
00199                 if ((strncmp(s, "ENDMDL", 6) == 0) && (ter_stop==0)){
00200                         ter_stop =1;
00201                         count_stop = count;
00202                 }
00203                 if (strncmp(s, "END", 6) == 0){
00204                         break;
00205                 }
00206                 //if ((strncmp(s, "TER", 3) == 0) && (ter_stop ==0)){
00207                         //if (count !=0) {ter_stop = 1;}
00208                         //count_stop = count;
00209                 //}
00210                 if (strncmp(s, "ATOM", 4) != 0) {
00211                         if (count == 0) {head.push_back(s);}
00212                         else {tail.push_back(s);}
00213                 }
00214                 else {
00215                         pWords.push_back(s);
00216                         atomName.push_back(pWords[count].substr(12,4));
00217                         residueName.push_back(pWords[count].substr(17,3));
00218                         chainId.push_back(pWords[count].substr(21,1));
00219                         elementSym.push_back(pWords[count].substr(76,2));
00220 
00221                         float x, y, z, tf;
00222                         int an, sn;
00223                         sscanf(&s[6], "%d", &an);
00224                         sscanf(&s[23], "%d %f %f %f %*f %f", &sn, &x, &y, &z, &tf);
00225 
00226                         if (count + 1 > get_number_points()) {
00227                                 set_number_points(2 * (count + 1)); 
00228                         }   //makes sure point array is big enough
00229 
00230                         points[4 * count] = x;
00231                         points[4 * count + 1] = y;
00232                         points[4 * count + 2] = z;
00233                         points[4 * count + 3] = tf;
00234                         pointInfo.push_back(an);
00235                         pointInfo.push_back(sn);
00236                         count++;
00237                 }
00238         }
00239 
00240         fclose(fp);
00241         set_number_points(count);
00242         return true;
00243 }
00244 
00245 
00246 //Accurately writes a pdb file once the information is alread stored.
00247 /*
00248 void PDBReader::save_to_pdb(const char *file) const {
00249         FILE *fp = fopen(file, "w");
00250         for (int i = 0; i <head.size(); i++) {
00251                 char head2 [200];
00252                 strcpy(head2, head[i].c_str());
00253                 fprintf (fp, head2);
00254         }
00255         int m = 0;
00256         for ( size_t i = 0; i < get_number_points(); i++) {
00257                 string curr = pWords[m];
00258                 string mid, final;
00259                 mid = curr.substr(12, 10);
00260                 final = curr.substr(76,2);
00261                 char mid2 [12];
00262                 strcpy(mid2, mid.c_str());
00263                 char final2 [4];
00264                 strcpy(final2, final.c_str());
00265                 m++;
00266                 fprintf(fp, "ATOM  %5d %10s%4d    %8.3f%8.3f%8.3f  1.00%6.2f          %2s\n", pointInfo[2*i], mid2, pointInfo[2*i +1], points[4 * i], points[4 * i + 1], points[4 * i + 2], points[4*i + 3], final2);
00267         }
00268         for (int i = 0; i <tail.size(); i++) {
00269                 char tail2 [200];
00270                 strcpy(tail2, tail[i].c_str());
00271                 fprintf (fp, tail2);
00272         }
00273         fprintf (fp, "END");
00274         fclose(fp);
00275 
00276 }
00277 */
00278 
00279 
00280 //Accurately writes a pdb file once the information is alread stored.
00281 void PDBReader::save_to_pdb(const char *file) const {
00282         FILE *fp = fopen(file, "w");
00283         int m = 0;
00284         for (size_t i =0; i< lines.size(); i++) {
00285                 char liner [200];
00286                 strcpy (liner, lines[i].c_str());
00287                 if (strncmp(liner, "ATOM", 4) != 0) {
00288                         fprintf (fp, "%s", liner);
00289                 }
00290                 else {
00291                         string curr = pWords[m];
00292                         string mid, final;
00293                         mid = curr.substr(12, 10);
00294                         final = curr.substr(76,2);
00295                         char mid2 [12];
00296                         strcpy(mid2, mid.c_str());
00297                         char final2 [4];
00298                         strcpy(final2, final.c_str());
00299                         fprintf(fp, "ATOM  %5d %10s%4d    %8.3f%8.3f%8.3f  1.00%6.2f          %2s\n", pointInfo[2*m], mid2, pointInfo[2*m +1], points[4 * m], points[4 * m + 1], points[4 * m + 2], points[4*m + 3], final2);
00300                         m++;
00301                 }
00302         }       
00303         fclose(fp);
00304 
00305 }
00306 
00307 
00308 //returns a vector of the x,y,z values 
00309 vector<float> PDBReader::get_points() {
00310 vector<float> ret;
00311 for (unsigned int i=0; i<n; i++) {
00312         ret.push_back((float)points[i*4]);
00313         ret.push_back((float)points[i*4+1]);
00314         ret.push_back((float)points[i*4+2]);
00315 }
00316 
00317 return ret;
00318 }
00319 
00320 //performs a right transform of all x,y,z points given a Transform object
00321 void PDBReader::right_transform(const Transform& transform) {
00322         for ( unsigned int i = 0; i < 4 * n; i += 4) {
00323                 Transform s = transform.transpose();
00324                 Vec3f v((float)points[i],(float)points[i+1],(float)points[i+2]);
00325                 v= s*v;
00326                 points[i]  =v[0];
00327                 points[i+1]=v[1];
00328                 points[i+2]=v[2];
00329         }
00330 }
00331 
00332 PointArray* PDBReader::makePointArray (const PDBReader& p) {
00333         PointArray* pArray  = new PointArray;
00334         p.save_to_pdb("thisFile3.txt");
00335         pArray->read_from_pdb("thisFile3.txt");
00336         remove("thisFile3.txt");
00337         
00338         return pArray;
00339 }
00340 
00341