Changeset 3846
- Timestamp:
- Jul 30, 2013 4:47:04 PM (11 years ago)
- Location:
- trunk/gui
- Files:
-
- 2 edited
-
scripts/drawing.tcl (modified) (1 diff)
-
src/RpPdbToVtk.c (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/gui/scripts/drawing.tcl
r3782 r3846 102 102 set pdbdata [$_xmlobj get $path.$elem.pdb] 103 103 if { $pdbdata != "" } { 104 if 1{104 if 0 { 105 105 global count 106 106 set f [open /tmp/file$count.pdb "w"] -
trunk/gui/src/RpPdbToVtk.c
r3782 r3846 21 21 #include "tcl.h" 22 22 23 24 static INLINE char * 25 SkipSpaces(char *string) 23 typedef struct { 24 int from, to; 25 } ConnectKey; 26 27 typedef struct { 28 int lineNum; 29 char atomName[6]; /* Atom name. */ 30 char symbolName[3]; /* Symbol name. */ 31 char residueName[4]; /* Residue name. */ 32 int serial; /* Serial # in PDB file. */ 33 double x, y, z; /* Coordinates of atom. */ 34 int number; /* Atomic number */ 35 } PdbAtom; 36 37 static int lineNum = 0; 38 static const char * 39 Itoa(int number) 40 { 41 static char buf[200]; 42 sprintf(buf, "%d", number); 43 return buf; 44 } 45 46 static INLINE const char * 47 SkipSpaces(const char *string) 26 48 { 27 49 while (isspace(*string)) { … … 31 53 } 32 54 33 static INLINE char * 34 GetLine(char **stringPtr, const char *endPtr, int *lengthPtr) 35 { 36 char *line, *p; 55 static int 56 IsSpaces(const char *string) 57 { 58 const char *p; 59 for (p = string; *p != '\0'; p++) { 60 if (!isspace(*p)) { 61 return 0; 62 } 63 } 64 return 1; 65 } 66 67 static INLINE const char * 68 GetLine(const char **stringPtr, const char *endPtr, int *lengthPtr) 69 { 70 const char *line; 71 const char *p; 37 72 38 73 line = SkipSpaces(*stringPtr); … … 47 82 return line; 48 83 } 49 50 static const char *elemNames[] = {51 "Hydrogen", /* 1 */52 "Helium", /* 2 */53 "Lithium", /* 3 */54 "Beryllium", /* 4 */55 "Boron", /* 5 */56 "Carbon", /* 6 */57 "Nitrogen", /* 7 */58 "Oxygen", /* 8 */59 "Fluorine", /* 9 */60 "Neon", /* 10 */61 "Sodium", /* 11 */62 "Magnesium", /* 12 */63 "Aluminium", /* 13 */64 "Silicon", /* 14 */65 "Phosphorus", /* 15 */66 "Sulfur", /* 16 */67 "Chlorine", /* 17 */68 "Argon", /* 18 */69 "Potassium", /* 19 */70 "Calcium", /* 20 */71 "Scandium", /* 21 */72 "Titanium", /* 22 */73 "Vanadium", /* 23 */74 "Chromium", /* 24 */75 "Manganese", /* 25 */76 "Iron", /* 26 */77 "Cobalt", /* 27 */78 "Nickel", /* 28 */79 "Copper", /* 29 */80 "Zinc", /* 30 */81 "Gallium", /* 31 */82 "Germanium", /* 32 */83 "Arsenic", /* 33 */84 "Selenium", /* 34 */85 "Bromine", /* 35 */86 "Krypton", /* 36 */87 "Rubidium", /* 37 */88 "Strontium", /* 38 */89 "Yttrium", /* 39 */90 "Zirconium", /* 40 */91 "Niobium", /* 41 */92 "Molybdenum", /* 42 */93 "Technetium", /* 43 */94 "Ruthenium", /* 44 */95 "Rhodium", /* 45 */96 "Palladium", /* 46 */97 "Silver", /* 47 */98 "Cadmium", /* 48 */99 "Indium", /* 49 */100 "Tin", /* 50 */101 "Antimony", /* 51 */102 "Tellurium", /* 52 */103 "Iodine", /* 53 */104 "Xenon", /* 54 */105 "Cesium", /* 55 */106 "Barium", /* 56 */107 "Lanthanum", /* 57 */108 "Cerium", /* 58 */109 "Praseodymium", /* 59 */110 "Neodymium", /* 60 */111 "Promethium", /* 61 */112 "Samarium", /* 62 */113 "Europium", /* 63 */114 "Gadolinium", /* 64 */115 "Terbium", /* 65 */116 "Dysprosium", /* 66 */117 "Holmium", /* 67 */118 "Erbium", /* 68 */119 "Thulium", /* 69 */120 "Ytterbium", /* 70 */121 "Lutetium", /* 71 */122 "Hafnium", /* 72 */123 "Tantalum", /* 73 */124 "Tungsten", /* 74 */125 "Rhenium", /* 75 */126 "Osmium", /* 76 */127 "Iridium", /* 77 */128 "Platinum", /* 78 */129 "Gold", /* 79 */130 "Mercury", /* 80 */131 "Thallium", /* 81 */132 "Lead", /* 82 */133 "Bismuth", /* 83 */134 "Polonium", /* 84 */135 "Astatine", /* 85 */136 "Radon", /* 86 */137 "Francium", /* 87 */138 "Radium", /* 88 */139 "Actinium", /* 89 */140 "Thorium", /* 90 */141 "Protactinium", /* 91 */142 "Uranium", /* 92 */143 "Neptunium", /* 93 */144 "Plutonium", /* 94 */145 "Americium", /* 95 */146 "Curium", /* 96 */147 "Berkelium", /* 97 */148 "Californium", /* 98 */149 "Einsteinium", /* 99 */150 "Fermium", /* 100 */151 "Mendelevium", /* 101 */152 "Nobelium", /* 102 */153 "Lawrencium", /* 103 */154 "Rutherfordium", /* 104 */155 "Dubnium", /* 105 */156 "Seaborgium", /* 106 */157 "Bohrium", /* 107 */158 "Hassium", /* 108 */159 "Meitnerium", /* 109 */160 "Darmstadtium", /* 110 */161 "Roentgenium", /* 111 */162 "Copernicium", /* 112 */163 "Ununtrium", /* 113 */164 "Flerovium", /* 114 */165 "Ununpentium", /* 115 */166 "Livermorium", /* 116 */167 "Ununseptium", /* 117 */168 "Ununoctium", /* 118 */169 NULL170 };171 84 172 85 static const char *symbolNames[] = { … … 293 206 294 207 int 295 VerifyElement(const char *elemName, const char *symbolName)208 GetAtomicNumber(Tcl_Interp *interp, PdbAtom *atomPtr) 296 209 { 297 210 const char **p; 211 char name[3], *namePtr; 298 212 int elemIndex, symbolIndex; 299 213 300 214 symbolIndex = elemIndex = -1; 301 elemName = SkipSpaces(elemName); 302 for (p = symbolNames; *p != NULL; p++) { 303 if (strcasecmp(elemName, *p) == 0) { 304 elemIndex = (p - symbolNames) + 1; 305 break; 306 } 307 } 308 symbolName = SkipSpaces(symbolName); 309 for (p = symbolNames; *p != NULL; p++) { 310 if (strcasecmp(symbolName, *p) == 0) { 311 symbolIndex = (p - symbolNames) + 1; 312 break; 215 216 /* 217 * The atomic symbol may be found the atom name in various locations 218 * "C " 219 * " C " 220 * " C " 221 " " C" 222 " "C 12" 223 " " C12" 224 */ 225 if ((atomPtr->atomName[0] == ' ') && (atomPtr->atomName[1] == ' ')) { 226 name[0] = atomPtr->atomName[2]; 227 name[1] = atomPtr->atomName[3]; 228 } else { 229 name[0] = atomPtr->atomName[0]; 230 name[1] = atomPtr->atomName[1]; 231 } 232 name[2] = '\0'; 233 if ((name[0] != ' ') || (name[1] != ' ')) { 234 namePtr = name; 235 if (name[0] == ' ') { 236 namePtr++; 237 } 238 if (name[1] == ' ') { 239 name[1] = '\0'; 240 } 241 for (p = symbolNames; *p != NULL; p++) { 242 if (strcasecmp(namePtr, *p) == 0) { 243 elemIndex = (p - symbolNames) + 1; 244 break; 245 } 246 } 247 } 248 name[0] = atomPtr->symbolName[0]; 249 name[1] = atomPtr->symbolName[1]; 250 name[2] = '\0'; 251 if (isdigit(name[1])) { 252 sscanf(name, "%d", &atomPtr->number); 253 return TCL_OK; 254 } 255 if ((name[0] != ' ') || (name[1] != ' ')) { 256 namePtr = name; 257 if (name[0] == ' ') { 258 namePtr++; 259 } 260 if (name[1] == ' ') { 261 name[1] = '\0'; 262 } 263 for (p = symbolNames; *p != NULL; p++) { 264 if (strcasecmp(namePtr, *p) == 0) { 265 symbolIndex = (p - symbolNames) + 1; 266 break; 267 } 313 268 } 314 269 } 315 270 if (symbolIndex > 0) { 316 271 if (elemIndex < 0) { 317 return symbolIndex; 272 atomPtr->number = symbolIndex; 273 return TCL_OK; 318 274 } 319 275 if (symbolIndex != elemIndex) { 320 return -1; 321 } 322 return elemIndex; 276 fprintf(stderr, "atomName %s and symbolName %s don't match\n", 277 atomPtr->atomName, atomPtr->symbolName); 278 atomPtr->number = symbolIndex; 279 return TCL_OK; 280 } 281 atomPtr->number = elemIndex; 282 return TCL_OK; 323 283 } else if (elemIndex > 0) { 324 return elemIndex; 325 } 326 return -1; 284 atomPtr->number = elemIndex; 285 return TCL_OK; 286 } 287 Tcl_AppendResult(interp, "line ", Itoa(lineNum), 288 ": bad atom \"", atomPtr->atomName, "\" or element \"", 289 atomPtr->symbolName, "\"", (char *)NULL); 290 return TCL_ERROR; 291 } 292 293 static int 294 SerialToIndex(Tcl_Interp *interp, Tcl_HashTable *tablePtr, const char *string, 295 int *indexPtr) 296 { 297 int serial; 298 long lserial, lindex; 299 Tcl_HashEntry *hPtr; 300 301 string = SkipSpaces(string); 302 if (Tcl_GetInt(interp, string, &serial) != TCL_OK) { 303 Tcl_AppendResult(interp, ": line ", Itoa(lineNum), 304 ": invalid serial number \"", string, 305 "\" in CONECT record", (char *)NULL); 306 return TCL_ERROR; 307 } 308 lserial = (long)serial; 309 hPtr = Tcl_FindHashEntry(tablePtr, (char *)lserial); 310 if (hPtr == NULL) { 311 Tcl_AppendResult(interp, "serial number \"", string, 312 "\" not found in table", (char *)NULL); 313 return TCL_ERROR; 314 } 315 lindex = (long)Tcl_GetHashValue(hPtr); 316 *indexPtr = (int)lindex; 317 return TCL_OK; 318 } 319 320 static int 321 GetCoordinates(Tcl_Interp *interp, const char *line, PdbAtom *atomPtr) 322 { 323 char buf[200]; 324 325 strncpy(buf, line + 30, 8); 326 buf[8] = '\0'; 327 if (Tcl_GetDouble(interp, buf, &atomPtr->x) != TCL_OK) { 328 Tcl_AppendResult(interp, "bad x-coordinate \"", buf, 329 "\"", (char *)NULL); 330 return TCL_ERROR; 331 } 332 strncpy(buf, line + 38, 8); 333 buf[8] = '\0'; 334 if (Tcl_GetDouble(interp, buf, &atomPtr->y) != TCL_OK) { 335 Tcl_AppendResult(interp, "bad y-coordinate \"", buf, 336 "\"", (char *)NULL); 337 return TCL_ERROR; 338 } 339 strncpy(buf, line + 46, 8); 340 buf[8] = '\0'; 341 if (Tcl_GetDouble(interp, buf, &atomPtr->z) != TCL_OK) { 342 Tcl_AppendResult(interp, "bad z-coordinate \"", buf, 343 "\"", (char *)NULL); 344 return TCL_ERROR; 345 } 346 return TCL_OK; 327 347 } 328 348 … … 330 350 * PdbToVtk string 331 351 */ 332 333 352 static int 334 353 PdbToVtkCmd(ClientData clientData, Tcl_Interp *interp, int objc, … … 336 355 { 337 356 Tcl_Obj *objPtr, *pointsObjPtr, *atomsObjPtr, *verticesObjPtr; 338 c har *p, *pend;357 const char *p, *pend; 339 358 char *string; 340 359 char mesg[2000]; 341 int length, n umAtoms, numConnections;342 Tcl_HashTable serialTable ;360 int length, nextIndex; 361 Tcl_HashTable serialTable, conectTable; 343 362 int i; 344 363 Tcl_HashEntry *hPtr; … … 346 365 int isNew; 347 366 348 numConnections = numAtoms= 0;367 lineNum = nextIndex = 0; 349 368 if (objc != 2) { 350 369 Tcl_AppendResult(interp, "wrong # arguments: should be \"", … … 353 372 } 354 373 Tcl_InitHashTable(&serialTable, TCL_ONE_WORD_KEYS); 374 Tcl_InitHashTable(&conectTable, sizeof(ConnectKey) / sizeof(int)); 355 375 string = Tcl_GetStringFromObj(objv[1], &length); 356 376 pointsObjPtr = Tcl_NewStringObj("", -1); … … 362 382 objPtr = NULL; 363 383 for (p = string, pend = p + length; p < pend; /*empty*/) { 364 c har *line, *q;384 const char *line; 365 385 char c; 366 386 int lineLength; 367 387 388 lineLength = 0; /* Suppress compiler warning. */ 389 lineNum++; 368 390 line = GetLine(&p, pend, &lineLength); 369 391 if (line >= pend) { … … 378 400 ((c == 'H') && (strncmp(line, "HETATM", 6) == 0))) { 379 401 char buf[200]; 380 char atomName[6]; 381 char symbolName[3]; 382 int serial; 383 double x, y, z; 384 int atom; 402 long lserial; 403 PdbAtom atom; 385 404 386 405 if (lineLength < 47) { … … 391 410 strncpy(buf, line + 6, 5); 392 411 buf[5] = '\0'; 393 if (Tcl_GetInt(interp, buf, &serial) != TCL_OK) { 394 Tcl_AppendResult(interp, "bad serial number \"", buf, 395 "\"", (char *)NULL); 396 goto error; 397 } 398 hPtr = Tcl_CreateHashEntry(&serialTable, (char *)serial, &isNew); 412 if (Tcl_GetInt(interp, buf, &atom.serial) != TCL_OK) { 413 Tcl_AppendResult(interp, "line ", Itoa(lineNum), 414 ": bad serial number \"", buf, "\"", (char *)NULL); 415 goto error; 416 } 417 lserial = (long)atom.serial; 418 hPtr = Tcl_CreateHashEntry(&serialTable, (char *)lserial, &isNew); 399 419 if (!isNew) { 400 Tcl_AppendResult(interp, "serial number \"", buf, 401 "\" found more than once", (char *)NULL); 402 goto error; 403 } 404 Tcl_SetHashValue(hPtr, numAtoms); 420 Tcl_AppendResult(interp, "line ", Itoa(lineNum), 421 ": serial number \"", buf, "\" found more than once", 422 (char *)NULL); 423 goto error; 424 } 425 Tcl_SetHashValue(hPtr, (long)nextIndex); 426 strncpy(atom.atomName, line + 12, 4); 427 atom.atomName[4] = '\0'; 428 strncpy(atom.residueName, line + 17, 3); 429 atom.residueName[3] = '\0'; 430 431 if (GetCoordinates(interp, line, &atom) != TCL_OK) { 432 goto error; 433 } 434 atom.symbolName[0] = '\0'; 435 if (lineLength >= 77) { 436 atom.symbolName[0] = line[76]; 437 atom.symbolName[1] = line[77]; 438 atom.symbolName[2] = '\0'; 439 } 440 if (GetAtomicNumber(interp, &atom) != TCL_OK) { 441 goto error; 442 } 405 443 Tcl_ListObjAppendElement(interp, verticesObjPtr, 406 444 Tcl_NewIntObj(1)); 407 445 Tcl_ListObjAppendElement(interp, verticesObjPtr, 408 Tcl_NewIntObj(serial)); 409 410 strncpy(atomName, line + 12, 4); 411 atomName[4] = '\0'; 412 413 strncpy(buf, line + 30, 8); 414 buf[8] = '\0'; 415 if (Tcl_GetDouble(interp, buf, &x) != TCL_OK) { 416 Tcl_AppendResult(interp, "bad x-coordinate \"", buf, 417 "\"", (char *)NULL); 418 goto error; 419 } 420 Tcl_ListObjAppendElement(interp, pointsObjPtr, Tcl_NewDoubleObj(x)); 421 strncpy(buf, line + 38, 8); 422 buf[8] = '\0'; 423 if (Tcl_GetDouble(interp, buf, &y) != TCL_OK) { 424 Tcl_AppendResult(interp, "bad y-coordinate \"", buf, 425 "\"", (char *)NULL); 426 goto error; 427 } 428 Tcl_ListObjAppendElement(interp, pointsObjPtr, Tcl_NewDoubleObj(y)); 429 strncpy(buf, line + 46, 8); 430 buf[8] = '\0'; 431 if (Tcl_GetDouble(interp, buf, &z) != TCL_OK) { 432 Tcl_AppendResult(interp, "bad z-coordinate \"", buf, 433 "\"", (char *)NULL); 434 goto error; 435 } 436 Tcl_ListObjAppendElement(interp, pointsObjPtr, Tcl_NewDoubleObj(z)); 437 symbolName[0] = '\0'; 438 if (lineLength >= 78) { 439 symbolName[0] = line[76]; 440 symbolName[1] = line[77]; 441 symbolName[2] = '\0'; 442 } 443 atom = VerifyElement(SkipSpaces(atomName), SkipSpaces(symbolName)); 444 if (atom < 0) { 445 Tcl_AppendResult(interp, "bad atom \"", atomName, 446 "\" or element \"", symbolName, "\"", (char *)NULL); 447 goto error; 448 } 449 Tcl_ListObjAppendElement(interp, atomsObjPtr, Tcl_NewIntObj(atom)); 450 numAtoms++; 446 Tcl_NewIntObj(atom.serial)); 447 Tcl_ListObjAppendElement(interp, pointsObjPtr, 448 Tcl_NewDoubleObj(atom.x)); 449 Tcl_ListObjAppendElement(interp, pointsObjPtr, 450 Tcl_NewDoubleObj(atom.y)); 451 Tcl_ListObjAppendElement(interp, pointsObjPtr, 452 Tcl_NewDoubleObj(atom.z)); 453 Tcl_ListObjAppendElement(interp, atomsObjPtr, 454 Tcl_NewIntObj(atom.number)); 455 nextIndex++; 451 456 } else if ((c == 'C') && (strncmp(line, "CONECT", 6) == 0)) { 452 numConnections++; 457 int a, i, n; 458 char buf[200]; 459 460 strncpy(buf, line + 6, 5); 461 buf[5] = '\0'; 462 if (lineLength < 11) { 463 char *msg; 464 465 msg = (char *)line; 466 msg[lineLength] = '\0'; 467 Tcl_AppendResult(interp, "line ", Itoa(lineNum), 468 ": bad CONECT record \"", msg, "\"", 469 (char *)NULL); 470 goto error; 471 } 472 if (SerialToIndex(interp, &serialTable, buf, &a) != TCL_OK) { 473 goto error; 474 } 475 /* Allow basic maximum of 4 connections. */ 476 for (n = 11, i = 0; i < 4; i++, n += 5) { 477 ConnectKey key; 478 int b; 479 if (n >= lineLength) { 480 break; 481 } 482 strncpy(buf, line + n, 5); 483 buf[5] = '\0'; 484 if (IsSpaces(buf)) { 485 break; /* No more entries */ 486 } 487 if (SerialToIndex(interp, &serialTable, buf, &b) != TCL_OK) { 488 goto error; 489 } 490 if (a > b) { 491 key.from = b; 492 key.to = a; 493 } else { 494 key.from = a; 495 key.to = b; 496 } 497 Tcl_CreateHashEntry(&conectTable, (char *)&key, &isNew); 498 } 453 499 } 454 500 } 455 501 objPtr = Tcl_NewStringObj("# vtk DataFile Version 2.0\n", -1); 456 Tcl_AppendToObj(objPtr, "Converted from PDB f ile\n", -1);502 Tcl_AppendToObj(objPtr, "Converted from PDB format\n", -1); 457 503 Tcl_AppendToObj(objPtr, "ASCII\n", -1); 458 504 Tcl_AppendToObj(objPtr, "DATASET POLYDATA\n", -1); 459 sprintf(mesg, "POINTS %d float\n", numAtoms);505 sprintf(mesg, "POINTS %d float\n", serialTable.numEntries); 460 506 Tcl_AppendToObj(objPtr, mesg, -1); 461 507 Tcl_AppendObjToObj(objPtr, pointsObjPtr); … … 463 509 Tcl_AppendToObj(objPtr, mesg, -1); 464 510 465 #ifdef notdef 466 if (numConnections > 0) { 467 sprintf(mesg, "LINES %d %d\n", numConnections, numConnections * 3);511 if (conectTable.numEntries > 0) { 512 sprintf(mesg, "LINES %d %d\n", conectTable.numEntries, 513 conectTable.numEntries * 3); 468 514 Tcl_AppendToObj(objPtr, mesg, -1); 469 for (hPtr = Tcl_FirstHashEntry(&con nectTable, &iter); hPtr != NULL;515 for (hPtr = Tcl_FirstHashEntry(&conectTable, &iter); hPtr != NULL; 470 516 hPtr = Tcl_NextHashEntry(&iter)) { 471 connectPtr = Tcl_GetHashEntry(hPtr); 472 sprintf(mesg, "2 %d %d\n", connectPtr->from, connectPtr->to); 473 } 474 } 475 Tcl_AppendObjToObj(objPtr, linesObjPtr); 476 sprintf(mesg, "\n"); 517 ConnectKey *keyPtr; 518 519 keyPtr = (ConnectKey *)Tcl_GetHashKey(&conectTable, hPtr); 520 sprintf(mesg, "2 %d %d\n", keyPtr->from, keyPtr->to); 521 Tcl_AppendToObj(objPtr, mesg, -1); 522 } 523 sprintf(mesg, "\n"); 524 Tcl_AppendToObj(objPtr, mesg, -1); 525 } 526 527 sprintf(mesg, "VERTICES %d %d\n", serialTable.numEntries, 528 serialTable.numEntries * 2); 477 529 Tcl_AppendToObj(objPtr, mesg, -1); 478 #endif 479 480 sprintf(mesg, "VERTICES %d %d\n", numAtoms, numAtoms * 2); 481 Tcl_AppendToObj(objPtr, mesg, -1); 482 for (i = 0; i < numAtoms; i++) { 530 for (i = 0; i < serialTable.numEntries; i++) { 483 531 sprintf(mesg, " 1 %d\n", i); 484 532 Tcl_AppendToObj(objPtr, mesg, -1); 485 533 } 486 sprintf(mesg, "POINT_DATA %d\n", numAtoms);534 sprintf(mesg, "POINT_DATA %d\n", serialTable.numEntries); 487 535 Tcl_AppendToObj(objPtr, mesg, -1); 488 536 sprintf(mesg, "SCALARS element int\n"); … … 492 540 Tcl_AppendObjToObj(objPtr, atomsObjPtr); 493 541 Tcl_DeleteHashTable(&serialTable); 542 Tcl_DeleteHashTable(&conectTable); 494 543 Tcl_DecrRefCount(pointsObjPtr); 495 544 Tcl_DecrRefCount(atomsObjPtr); … … 500 549 return TCL_OK; 501 550 error: 551 Tcl_DeleteHashTable(&conectTable); 502 552 Tcl_DeleteHashTable(&serialTable); 503 553 Tcl_DecrRefCount(pointsObjPtr);
Note: See TracChangeset
for help on using the changeset viewer.
