/* * pdb.c * * Summary: Read trajectories from PDB file. * * Author: David Hardy */ #include #include #include #include "mdtypes.h" #include "pdb.h" #include "pdbrec.h" #include "error.h" int pdb_startup(void) { return pdbrec_startup(); } void pdb_shutdown(void) { pdbrec_shutdown(); } int pdb_init(Pdb *p) { memset(p, 0, sizeof(Pdb)); return 0; } void pdb_destroy(Pdb *p) { /* nothing to do! */ } int pdb_open(Pdb *p, const char *filename, int filetype) { p->f = fopen(filename, "r"); MDIO_ASSERT(p->f != NULL); return pdbrec_init(&(p->pdbrec), filetype); } int pdb_close(Pdb *p) { int ret = fclose(p->f); MDIO_ASSERT(ret == 0); memset(p, 0, sizeof(Pdb)); return ret; } int pdb_read(Pdb *p, void *vec, void *a, void *b, int elemtype, int buflen, int *numread) { int mask, i, k, ret; FILE *f = p->f; Pdbrec *pdbrec = &(p->pdbrec); Pdbrec_Data *data = &(p->data); MD_Fvec *fvec = NULL; float *fa = NULL, *fb = NULL; MD_Dvec *dvec = NULL; double *da = NULL, *db = NULL; char buf[128]; mask = -1; /* block everything */ if (vec != NULL) { UNMASK(mask, 8); UNMASK(mask, 9); UNMASK(mask, 10); } if (a != NULL) { UNMASK(mask, 11); } if (b != NULL) { UNMASK(mask, 12); } p->data.d.atom.mask = mask; MDIO_ASSERT(elemtype == FLOATTYPE || elemtype == DOUBLETYPE); if (elemtype == FLOATTYPE) { fvec = (MD_Fvec *) vec; fa = (float *) a; fb = (float *) b; } else { dvec = (MD_Dvec *) vec; da = (double *) a; db = (double *) b; } k = 0; while (k < buflen && fgets(buf, sizeof(buf), f) != NULL) { /* pad out string with blanks */ i = strlen(buf); if (buf[i-1] == '\n') i--; while (i < 80) buf[i++] = ' '; buf[80] = '\n'; buf[81] = '\0'; MDIO_ASSERT(strlen(buf) == 81); /* 80 chars + newline */ ret = pdbrec_extract(pdbrec, buf); MDIO_ASSERT(ret >= 0 && ret <= PDBREC_END); if (ret == PDBREC_ATOM || ret == PDBREC_HETATM) { ret = pdbrec_getdata(pdbrec, data); MDIO_ASSERT(ret == 0); if (dvec != NULL) { dvec[k].x = (double) data->d.atom.x; dvec[k].y = (double) data->d.atom.y; dvec[k].z = (double) data->d.atom.z; } else if (fvec != NULL) { fvec[k].x = data->d.atom.x; fvec[k].y = data->d.atom.y; fvec[k].z = data->d.atom.z; } if (da != NULL) { da[k] = (double) data->d.atom.occupancy; } else if (fa != NULL) { fa[k] = data->d.atom.occupancy; } if (db != NULL) { db[k] = (double) data->d.atom.tempFactor; } else if (fb != NULL) { fb[k] = data->d.atom.tempFactor; } k++; } } *numread = k; return 0; }