00001
00007
00008
00009
00010
00011
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
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
00042 while ( fgets(buf, 150, infile) ) {
00043 PDBData *newelement;
00044 char *s;
00045 for (s=buf; *s && *s!='\n'; s++)
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;
00054
00055
00056 newelement = new_PDBData(buf);
00057 if (!newelement) {
00058 NAMD_die("Could not allocate PDBData.\n");
00059 }
00060
00061
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 }
00078 Fclose(infile);
00079 }
00080 #endif
00081
00082
00083 PDB::PDB( const char *pdbfilename) {
00084 FILE *infile;
00085 char buf[160];
00086
00087 atomCount = 0;
00088 atomListHead = atomListTail = NULL;
00089 infile = Fopen(pdbfilename, "r");
00090 if (! infile) {
00091 char s[500];
00092 sprintf(s, "Cannot open file '%s' for input in PDB::PDB.", pdbfilename);
00093 NAMD_err(s);
00094 }
00095
00096
00097 while ( fgets(buf, 150, infile) ) {
00098 PDBData *newelement;
00099 char *s;
00100 for (s=buf; *s && *s!='\n'; s++)
00101 ;
00102 *s = 0;
00103 if ( s == (buf + 149) ) {
00104 char s[500];
00105 sprintf( s, "Line too long in pdbfile %s:\n%s\n", pdbfilename, buf);
00106 NAMD_die(s);
00107 }
00108 *(s+1) = 0;
00109
00110
00111 newelement = new_PDBData(buf);
00112 if (!newelement) {
00113 NAMD_die("Could not allocate PDBData.\n");
00114 }
00115
00116
00117 if (newelement -> type() != PDBData::ATOM &&
00118 newelement -> type() != PDBData::HETATM) {
00119 delete newelement;
00120 } else {
00121 add_atom_element( (PDBAtom *) newelement);
00122 }
00123 }
00124 Fclose(infile);
00125
00126
00127
00128
00129 {
00130
00131 #ifdef MEM_OPT_VERSION
00132
00133 atomArray = new PDBCoreData[atomCount];
00134 if(atomArray == NULL)
00135 NAMD_die("memory allocation failed in PDB::PDB");
00136
00137 PDBAtomList *tmp = atomListHead;
00138 for(int i=0; tmp!=NULL; tmp=tmp->next, i++){
00139 const BigReal *coor = tmp->data->coordinates();
00140 atomArray[i].coor[0] = coor[0];
00141 atomArray[i].coor[1] = coor[1];
00142 atomArray[i].coor[2] = coor[2];
00143 atomArray[i].myoccupancy = tmp->data->occupancy();
00144 atomArray[i].tempfactor = tmp->data->temperaturefactor();
00145 }
00146
00147
00148 tmp = atomListHead;
00149 while(tmp!=NULL){
00150 PDBAtomList *nextTmp = tmp->next;
00151 delete tmp->data;
00152 delete tmp;
00153 tmp = nextTmp;
00154 }
00155 atomListHead = atomListTail = NULL;
00156 #else
00157 atomArray = new PDBAtomPtr[atomCount];
00158 if ( atomArray == NULL )
00159 {
00160 NAMD_die("memory allocation failed in PDB::PDB");
00161 }
00162 PDBAtomList *tmp = atomListHead;
00163 int i=0;
00164 for (i=0, tmp = atomListHead; tmp != NULL; tmp = tmp -> next, i++) {
00165 atomArray[i] = tmp -> data;
00166 }
00167
00168 PDBAtomList *tmp2;
00169 for (tmp2 = tmp = atomListHead; tmp != NULL; tmp = tmp2) {
00170 tmp2 = tmp->next;
00171 delete tmp;
00172 }
00173 atomListHead = atomListTail = NULL;
00174 #endif
00175 }
00176
00177 }
00178
00179
00180
00181 PDB::~PDB( void )
00182 {
00183 #ifndef MEM_OPT_VERSION
00184 int i;
00185 for (i=atomCount-1; i>=0; i--)
00186 delete atomArray[i];
00187 #endif
00188 delete [] atomArray;
00189 atomArray = NULL;
00190 atomCount = 0;
00191 }
00192
00193
00194 void PDB::write(const char *outfilename, const char *commentline)
00195 {
00196 int i;
00197 char s[200];
00198 FILE *outfile;
00199 if ((outfile = fopen(outfilename, "w")) == NULL) {
00200 sprintf(s, "Cannot open file '%s' in PDB::write.", outfilename);
00201 NAMD_err(s);
00202 }
00203
00204 if (commentline != NULL)
00205 {
00206 sprintf(s, "REMARK %s\n", commentline);
00207 if (fputs(s, outfile) == EOF)
00208 {
00209 NAMD_err("EOF in PDB::write writing the comment line - file system full?");
00210 }
00211 }
00212
00213 for (i=0; i<atomCount; i++){
00214 #ifdef MEM_OPT_VERSION
00215 atomArray[i].sprint(s, PDBData::COLUMNS);
00216 #else
00217 atomArray[i]->sprint(s, PDBData::COLUMNS);
00218 #endif
00219 if ( (fputs(s, outfile) == EOF) ||
00220 (fputs("\n", outfile) == EOF) ) {
00221 sprintf(s, "EOF in PDB::write line %d - file system full?", i);
00222 NAMD_err(s);
00223 }
00224 }
00225 if (fputs("END\n", outfile) == EOF) {
00226 NAMD_err("EOF in PDB::write while printing 'END' -- file system full?");
00227 }
00228 if (fclose(outfile) == EOF) {
00229 NAMD_err("EOF in PDB::write while closing -- file system full?");
00230 }
00231
00232 }
00233
00234
00235 void PDB::add_atom_element( PDBAtom *newAtom)
00236 {
00237 PDBAtomList *tmp = new PDBAtomList;
00238 if ( tmp == NULL )
00239 {
00240 NAMD_die("memory allocation failed in PDB::add_atom_element");
00241 }
00242 tmp -> data = newAtom;
00243 tmp -> next = NULL;
00244
00245 if (atomListHead == NULL) {
00246 atomListHead = atomListTail = tmp;
00247 } else {
00248 atomListTail -> next = tmp;
00249 atomListTail = tmp;
00250 }
00251 atomCount++;
00252 }
00253
00254
00255
00256 int PDB::num_atoms( void)
00257 {
00258 return atomCount;
00259 }
00260
00261
00262
00263
00264 void PDB::set_all_positions(Vector *pos)
00265 {
00266 #ifdef MEM_OPT_VERSION
00267 for(int i=0; i<atomCount; i++){
00268 atomArray[i].coor[0] = pos[i].x;
00269 atomArray[i].coor[1] = pos[i].y;
00270 atomArray[i].coor[2] = pos[i].z;
00271 }
00272 #else
00273 int i;
00274 PDBAtomPtr *atomptr;
00275
00276 for (i=0, atomptr=atomArray; i<atomCount; atomptr++, i++)
00277 {
00278 (*atomptr)->xcoor(pos[i].x);
00279 (*atomptr)->ycoor(pos[i].y);
00280 (*atomptr)->zcoor(pos[i].z);
00281 }
00282 #endif
00283 }
00284
00285 void PDB::get_position_for_atom(Vector *pos, int aid){
00286 #ifdef MEM_OPT_VERSION
00287 pos->x = atomArray[aid].coor[0];
00288 pos->y = atomArray[aid].coor[1];
00289 pos->z = atomArray[aid].coor[2];
00290 #else
00291 PDBAtomPtr *atomptr = &atomArray[aid];
00292 pos->x = (*atomptr)->xcoor();
00293 pos->y = (*atomptr)->ycoor();
00294 pos->z = (*atomptr)->zcoor();
00295 #endif
00296 }
00297
00298
00299 void PDB::get_all_positions(Vector *pos)
00300 {
00301 #ifdef MEM_OPT_VERSION
00302 for(int i=0; i<atomCount; i++){
00303 pos[i].x = atomArray[i].coor[0];
00304 pos[i].y = atomArray[i].coor[1];
00305 pos[i].z = atomArray[i].coor[2];
00306 }
00307 #else
00308 int i;
00309 PDBAtomPtr *atomptr;
00310
00311 for (i=0, atomptr=atomArray; i<atomCount; atomptr++, i++)
00312 {
00313 pos[i].x = (*atomptr)->xcoor();
00314 pos[i].y = (*atomptr)->ycoor();
00315 pos[i].z = (*atomptr)->zcoor();
00316 }
00317 #endif
00318 }
00319
00320
00321 #ifdef MEM_OPT_VERSION
00322 PDBCoreData *PDB::atom(int place){
00323 if(place<0 || place>=atomCount) return NULL;
00324 return &atomArray[place];
00325 }
00326 #else
00327 PDBAtom *PDB::atom(int place)
00328 {
00329 if (place <0 || place >= atomCount)
00330 return NULL;
00331 return atomArray[place];
00332 }
00333 #endif
00334
00335
00336
00337 void PDB::find_extremes(BigReal *min, BigReal *max, Vector rec, BigReal frac) const
00338 {
00339 SortableResizeArray<BigReal> coor;
00340 coor.resize(atomCount);
00341 SortableResizeArray<BigReal>::iterator c_i = coor.begin();
00342 #ifdef MEM_OPT_VERSION
00343 PDBCoreData *atomptr = atomArray;
00344 for(int i=0; i<atomCount; i++, atomptr++){
00345 Vector pos(atomptr->xcoor(),atomptr->ycoor(),atomptr->zcoor());
00346 c_i[i] = rec*pos;
00347 }
00348 #else
00349 PDBAtomPtr *atomptr = atomArray;
00350 for (int i=0; i<atomCount; ++i, ++atomptr) {
00351 PDBAtom *atom = *atomptr;
00352 Vector pos(atom->xcoor(),atom->ycoor(),atom->zcoor());
00353 c_i[i] = rec*pos;
00354 }
00355 #endif
00356 coor.sort();
00357 int ilow = (int)((1.0 - frac) * atomCount);
00358 if ( ilow < 0 ) ilow = 0;
00359 if ( ilow > atomCount/2 ) ilow = atomCount/2;
00360 int ihigh = atomCount - ilow - 1;
00361 BigReal span = coor[ihigh] - coor[ilow];
00362 BigReal extension = (1.0 - frac) * span / (2.0 * frac - 1.0);
00363 *max = coor[ihigh] + extension;
00364 *min = coor[ilow] - extension;
00365 }
00366
00367
00368 #ifdef TEST_PDB_CLASS
00369
00370 main()
00371 {
00372 PDB *pdb = new PDB("pti.pdb");
00373 if ( atomArray == NULL )
00374 {
00375 NAMD_die("memory allocation failed in main of test PDB class");
00376 }
00377 ilist = pdb->find_atom_name("CA");
00378 if (!ilist)
00379 printf("None found.\n");
00380 else {
00381 int i;
00382 char s[200];
00383 for (i=0; i<ilist->num(); i++) {
00384 pdb->atom((*ilist)[i]) -> sprint(s);
00385 printf("%s\n", s);
00386 }
00387 delete ilist;
00388 }
00389
00390 printf("Now, search through space.\n");
00391
00392 ilist = pdb->find_atoms_in_region(4.38, 19.5, 3.0,
00393 4.40, 20.0, 3.2 );
00394 if (!ilist)
00395 printf("None found.\n");
00396 else {
00397 int i;
00398 char s[200];
00399 printf("%d found\n", ilist -> num());
00400 for (i=0; i<ilist->num(); i++) {
00401 pdb->atom((*ilist)[i]) -> sprint(s);
00402 printf("%s\n", s);
00403 }
00404 delete ilist;
00405 }
00406 }
00407 #endif // TEST_PDB_CLASS
00408
00409
00410
00411
00412
00413 static int readtoeoln(FILE *f) {
00414 char c;
00415
00416
00417 while((c = getc(f)) != '\n') {
00418 if (c == EOF)
00419 return -1;
00420 }
00421
00422 return 0;
00423 }
00424
00425
00426
00427
00428 PDB::PDB( const char *filename, Ambertoppar *amber_data)
00429 { int i,j,k;
00430 BigReal coor[3];
00431 char buf[13],resname[5],atomname[5];
00432 FILE *infile;
00433 PDBAtom *pdb;
00434
00435 if ((infile=Fopen(filename, "r")) == NULL)
00436 NAMD_err("Can't open AMBER coordinate file!");
00437
00438 readtoeoln(infile);
00439
00440 fscanf(infile,"%d",&atomCount);
00441 if (atomCount != amber_data->Natom)
00442 NAMD_die("Num of atoms in coordinate file is different from that in parm file!");
00443 readtoeoln(infile);
00444
00445 #ifdef MEM_OPT_VERSION
00446 atomArray = new PDBCoreData[atomCount];
00447 #else
00448 atomArray = new PDBAtomPtr[atomCount];
00449 #endif
00450
00451 if ( atomArray == NULL )
00452 {
00453 NAMD_die("memory allocation failed in PDB::PDB");
00454 }
00455
00456
00457
00458 for (i=0; i<atomCount; ++i)
00459 {
00460 for (j=0; j<3; ++j)
00461 { for (k=0; k<12; ++k)
00462 { buf[k]=getc(infile);
00463 if (buf[k]=='\n' || buf[k]=='\0' || buf[k]==EOF)
00464 NAMD_die("Error reading AMBER coordinate file!");
00465 }
00466 buf[12] = '\0';
00467 coor[j] = atof(buf);
00468 }
00469 if (i%2 == 1)
00470 readtoeoln(infile);
00471 #ifdef MEM_OPT_VERSION
00472 atomArray[i].coor[0] = coor[0];
00473 atomArray[i].coor[1] = coor[1];
00474 atomArray[i].coor[2] = coor[2];
00475 atomArray[i].myoccupancy = PDBAtom::default_occupancy;
00476 atomArray[i].tempfactor = PDBAtom::default_temperaturefactor;
00477 #else
00478
00479 for (j=0; j<4; ++j)
00480 { resname[j] = amber_data->ResNames[amber_data->AtomRes[i]*4+j];
00481 atomname[j] = amber_data->AtomNames[i*4+j];
00482 }
00483 resname[4] = atomname[4] = '\0';
00484
00485 pdb = new PDBAtomRecord("");
00486 pdb->name(atomname);
00487 pdb->residuename(resname);
00488 pdb->serialnumber(i+1);
00489 pdb->residueseq(amber_data->AtomRes[i]+1);
00490 pdb->coordinates(coor);
00491 atomArray[i] = pdb;
00492 #endif
00493 }
00494 }
00495
00496 #define LINESIZE 100
00497
00498
00499
00500
00501
00502 PDB::PDB(const char *filename, const GromacsTopFile *topology) {
00503 int i;
00504 char buf[LINESIZE];
00505 FILE *infile;
00506
00507
00508 infile=Fopen(filename, "r");
00509 if (infile == NULL)
00510 NAMD_err("Can't open GROMACS coordinate file!");
00511
00512 fgets(buf,LINESIZE-1,infile);
00513
00514
00515
00516 fgets(buf,LINESIZE-1,infile);
00517 sscanf(buf,"%d",&atomCount);
00518 if (atomCount != topology->getNumAtoms())
00519 NAMD_die("Num of atoms in coordinate file is different from that in topology file!");
00520
00521
00522 #ifdef MEM_OPT_VERSION
00523 atomArray = new PDBCoreData[atomCount];
00524 #else
00525 atomArray = new PDBAtomPtr[atomCount];
00526 #endif
00527 if ( atomArray == NULL )
00528 NAMD_die("memory allocation failed in PDB::PDB");
00529
00530 #ifdef MEM_OPT_VERSION
00531 for(i=0; i<atomCount; i++){
00532 fgets(buf,LINESIZE-1,infile);
00533 char *buf2 = buf+20;
00534 BigReal coor[3];
00535 if(3 != sscanf(buf2,"%lf%lf%lf", &coor[0],&coor[1],&coor[2]))
00536 NAMD_die("Couldn't get three coordinates from file.");
00537
00538 coor[0] *= 10;
00539 coor[1] *= 10;
00540 coor[2] *= 10;
00541
00542 atomArray[i].coor[0] = coor[0];
00543 atomArray[i].coor[1] = coor[1];
00544 atomArray[i].coor[2] = coor[2];
00545 atomArray[i].myoccupancy = PDBAtom::default_occupancy;
00546 atomArray[i].tempfactor = PDBAtom::default_temperaturefactor;
00547 }
00548 #else
00549 for (i=0;i<atomCount;i++) {
00550 char *buf2, resname[11], atomname[11], atmtype[11];
00551 int resnum, typenum;
00552 Real charge,mass;
00553 BigReal coor[3];
00554 PDBAtom *pdb = new PDBAtomRecord("");
00555
00556 fgets(buf,LINESIZE-1,infile);
00557 buf2 = buf+20;
00558 if(3 != sscanf(buf2,"%lf%lf%lf",
00559 &coor[0],&coor[1],&coor[2]))
00560 NAMD_die("Couldn't get three coordinates from file.");
00561 topology->getAtom(i,&resnum,resname,
00562 atomname,atmtype,&typenum,&charge,&mass);
00563 coor[0] *= 10;
00564 coor[1] *= 10;
00565 coor[2] *= 10;
00566
00567 pdb->name(atomname);
00568 pdb->residuename(resname);
00569 pdb->serialnumber(i+1);
00570 pdb->residueseq(resnum+1);
00571 pdb->coordinates(coor);
00572
00573 atomArray[i] = pdb;
00574 }
00575 #endif
00576 }
00577
00578 #ifdef MEM_OPT_VERSION
00579 void PDBCoreData::sprint( char *outstr, PDBData::PDBFormatStyle usestyle){
00580 if(usestyle == PDBData::COLUMNS){
00581
00582 for(int i=0; i<79; i++) outstr[i] = 32;
00583 outstr[79] = 0;
00584 PDBData::sprintcol(outstr, PDBAtom::SX, PDBAtom::LCOOR, PDBAtom::LCOORPREC, xcoor());
00585 PDBData::sprintcol(outstr, PDBAtom::SY, PDBAtom::LCOOR, PDBAtom::LCOORPREC, ycoor());
00586 PDBData::sprintcol(outstr, PDBAtom::SZ, PDBAtom::LCOOR, PDBAtom::LCOORPREC, zcoor());
00587 PDBData::sprintcol(outstr, PDBAtom::SOCC, PDBAtom::LOCC, PDBAtom::LOCCPREC, occupancy());
00588 PDBData::sprintcol(outstr, PDBAtom::STEMPF, PDBAtom::LTEMPF, PDBAtom::LTEMPFPREC, temperaturefactor());
00589 }else{
00590
00591 char tmpstr[50];
00592 if (xcoor() == PDBAtom::default_coor)
00593 sprintf(tmpstr, " #");
00594 else
00595 sprintf(tmpstr, " %*.*f", PDBAtom::LCOOR, PDBAtom::LCOORPREC, xcoor());
00596 strcat(outstr, tmpstr);
00597 if (ycoor() == PDBAtom::default_coor)
00598 sprintf(tmpstr, " #");
00599 else
00600 sprintf(tmpstr, " %*.*f", PDBAtom::LCOOR, PDBAtom::LCOORPREC, ycoor());
00601 strcat(outstr, tmpstr);
00602 if (zcoor() == PDBAtom::default_coor)
00603 sprintf(tmpstr, " #");
00604 else
00605 sprintf(tmpstr, " %*.*f", PDBAtom::LCOOR, PDBAtom::LCOORPREC, zcoor());
00606 strcat(outstr, tmpstr);
00607
00608 sprintf(tmpstr, " %*.*f", PDBAtom::LOCC, PDBAtom::LOCCPREC, occupancy());
00609 strcat(outstr, tmpstr);
00610
00611 sprintf(tmpstr, " %*.*f", PDBAtom::LTEMPF, PDBAtom::LTEMPFPREC, temperaturefactor());
00612 strcat(outstr, tmpstr);
00613 }
00614 }
00615 #endif
00616