00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
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
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
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
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
00143 vector<float> PDBReader::get_y() {
00144
00145 return y;
00146 }
00147
00148
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
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));
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
00207
00208
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 }
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
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
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
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
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