source: trunk/gui/src/RpPdbToVtk.c @ 4979

Last change on this file since 4979 was 4658, checked in by ldelgass, 9 years ago

fix warning
File size: 20.8 KB
Line 
1
2/*
3 * ----------------------------------------------------------------------
4 *  RpPdbToVtk -
5 *
6 * ======================================================================
7 *  AUTHOR:  Michael McLennan, Purdue University
8 *  Copyright (c) 2004-2012  HUBzero Foundation, LLC
9 *
10 *  See the file "license.terms" for information on usage and
11 *  redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES.
12 * ======================================================================
13 */
14#include <string.h>
15#include <stdlib.h>
16#include <ctype.h>
17#include <math.h>
18#include <limits.h>
19#include <float.h>
20#include "tcl.h"
21
22#define BOND_NONE       (0)
23#define BOND_CONECT     (1<<0)
24#define BOND_COMPUTE    (1<<1)
25#define BOND_BOTH       (BOND_CONECT|BOND_COMPUTE)
26
27typedef struct {
28    int from, to;
29} ConnectKey;
30
31typedef struct {
32    double x, y, z;                     /* Coordinates of atom. */
33    int number;                         /* Atomic number */
34    int ordinal;                        /* Index of the atom in VTK. */
35    int numConnections;                 /*  */
36} PdbAtom;
37
38typedef struct {
39    const char *symbol;                 /* Atom symbol. */
40    float radius;                       /* Covalent radius of atom. */
41} PeriodicElement;
42
43static PeriodicElement elements[] = {
44    {    "X",       0.1f,       },      /* 0 */
45    {    "H",       0.32f,      },      /* 1 */         
46    {    "He",      0.93f,      },      /* 2 */         
47    {    "Li",      1.23f,      },      /* 3 */         
48    {    "Be",      0.90f,      },      /* 4 */         
49    {    "B",       0.82f,      },      /* 5 */         
50    {    "C",       0.77f,      },      /* 6 */         
51    {    "N",       0.75f,      },      /* 7 */         
52    {    "O",       0.73f,      },      /* 8 */         
53    {    "F",       0.72f,      },      /* 9 */         
54    {    "Ne",      0.71f,      },      /* 10 */       
55    {    "Na",      1.54f,      },      /* 11 */       
56    {    "Mg",      1.36f,      },      /* 12 */       
57    {    "Al",      1.18f,      },      /* 13 */       
58    {    "Si",      1.11f,      },      /* 14 */       
59    {    "P",       1.06f,      },      /* 15 */       
60    {    "S",       1.02f,      },      /* 16 */       
61    {    "Cl",      0.99f,      },      /* 17 */       
62    {    "Ar",      0.98f,      },      /* 18 */       
63    {    "K",       2.03f,      },      /* 19 */       
64    {    "Ca",      1.74f,      },      /* 20 */       
65    {    "Sc",      1.44f,      },      /* 21 */       
66    {    "Ti",      1.32f,      },      /* 22 */       
67    {    "V",       1.22f,      },      /* 23 */       
68    {    "Cr",      1.18f,      },      /* 24 */       
69    {    "Mn",      1.17f,      },      /* 25 */       
70    {    "Fe",      1.17f,      },      /* 26 */       
71    {    "Co",      1.16f,      },      /* 27 */       
72    {    "Ni",      1.15f,      },      /* 28 */       
73    {    "Cu",      1.17f,      },      /* 29 */       
74    {    "Zn",      1.25f,      },      /* 30 */       
75    {    "Ga",      1.26f,      },      /* 31 */       
76    {    "Ge",      1.22f,      },      /* 32 */       
77    {    "As",      1.20f,      },      /* 33 */       
78    {    "Se",      1.16f,      },      /* 34 */       
79    {    "Br",      1.14f,      },      /* 35 */       
80    {    "Kr",      1.12f,      },      /* 36 */       
81    {    "Rb",      2.16f,      },      /* 37 */       
82    {    "Sr",      1.91f,      },      /* 38 */       
83    {    "Y",       1.62f,      },      /* 39 */       
84    {    "Zr",      1.45f,      },      /* 40 */       
85    {    "Nb",      1.34f,      },      /* 41 */       
86    {    "Mo",      1.30f,      },      /* 42 */       
87    {    "Tc",      1.27f,      },      /* 43 */       
88    {    "Ru",      1.25f,      },      /* 44 */       
89    {    "Rh",      1.25f,      },      /* 45 */       
90    {    "Pd",      1.28f,      },      /* 46 */       
91    {    "Ag",      1.34f,      },      /* 47 */       
92    {    "Cd",      1.48f,      },      /* 48 */       
93    {    "In",      1.44f,      },      /* 49 */       
94    {    "Sn",      1.41f,      },      /* 50 */       
95    {    "Sb",      1.40f,      },      /* 51 */       
96    {    "Te",      1.36f,      },      /* 52 */       
97    {    "I",       1.33f,      },      /* 53 */       
98    {    "Xe",      1.31f,      },      /* 54 */       
99    {    "Cs",      2.35f,      },      /* 55 */       
100    {    "Ba",      1.98f,      },      /* 56 */       
101    {    "La",      1.69f,      },      /* 57 */       
102    {    "Ce",      1.65f,      },      /* 58 */       
103    {    "Pr",      1.65f,      },      /* 59 */       
104    {    "Nd",      1.64f,      },      /* 60 */       
105    {    "Pm",      1.63f,      },      /* 61 */       
106    {    "Sm",      1.62f,      },      /* 62 */       
107    {    "Eu",      1.85f,      },      /* 63 */       
108    {    "Gd",      1.61f,      },      /* 64 */       
109    {    "Tb",      1.59f,      },      /* 65 */       
110    {    "Dy",      1.59f,      },      /* 66 */       
111    {    "Ho",      1.58f,      },      /* 67 */       
112    {    "Er",      1.57f,      },      /* 68 */       
113    {    "Tm",      1.56f,      },      /* 69 */       
114    {    "Yb",      1.74f,      },      /* 70 */       
115    {    "Lu",      1.56f,      },      /* 71 */       
116    {    "Hf",      1.44f,      },      /* 72 */       
117    {    "Ta",      1.34f,      },      /* 73 */       
118    {    "W",       1.30f,      },      /* 74 */       
119    {    "Re",      1.28f,      },      /* 75 */       
120    {    "Os",      1.26f,      },      /* 76 */       
121    {    "Ir",      1.27f,      },      /* 77 */       
122    {    "Pt",      1.30f,      },      /* 78 */       
123    {    "Au",      1.34f,      },      /* 79 */       
124    {    "Hg",      1.49f,      },      /* 80 */       
125    {    "Tl",      1.48f,      },      /* 81 */       
126    {    "Pb",      1.47f,      },      /* 82 */       
127    {    "Bi",      1.46f,      },      /* 83 */       
128    {    "Po",      1.46f,      },      /* 84 */       
129    {    "At",      1.45f,      },      /* 85 */       
130    {    "Rn",      1.50f,      },      /* 86 */       
131    {    "Fr",      2.60f,      },      /* 87 */       
132    {    "Ra",      2.21f,      },      /* 88 */       
133    {    "Ac",      2.15f,      },      /* 89 */       
134    {    "Th",      2.06f,      },      /* 90 */       
135    {    "Pa",      2.00f,      },      /* 91 */       
136    {    "U",       1.96f,      },      /* 92 */       
137    {    "Np",      1.90f,      },      /* 93 */       
138    {    "Pu",      1.87f,      },      /* 94 */       
139    {    "Am",      1.80f,      },      /* 95 */       
140    {    "Cm",      1.69f,      },      /* 96 */       
141    {    "Bk",      0.1f,       },      /* 97 */       
142    {    "Cf",      0.1f,       },      /* 98 */       
143    {    "Es",      0.1f,       },      /* 99 */       
144    {    "Fm",      0.1f,       },      /* 100 */       
145    {    "Md",      0.1f,       },      /* 101 */       
146    {    "No",      0.1f,       },      /* 102 */       
147    {    "Lr",      0.1f,       },      /* 103 */       
148    {    "Rf",      0.1f,       },      /* 104 */       
149    {    "Db",      0.1f,       },      /* 105 */       
150    {    "Sg",      0.1f,       },      /* 106 */       
151    {    "Bh",      0.1f,       },      /* 107 */       
152    {    "Hs",      0.1f,       },      /* 108 */       
153    {    "Mt",      0.1f,       },      /* 109 */       
154    {    "Ds",      0.1f,       },      /* 110 */       
155    {    "Rg",      0.1f,       },      /* 111 */       
156    {    "Cn",      0.1f,       },      /* 112 */       
157    {    "Uut",     0.1f,       },      /* 113 */       
158    {    "Fl",      0.1f,       },      /* 114 */       
159    {    "Uup",     0.1f,       },      /* 115 */       
160    {    "Lv",      0.1f,       },      /* 116 */       
161    {    "Uus",     0.1f,       },      /* 117 */       
162    {    "Uuo",     0.1f,       }       /* 118 */       
163};
164
165static int numPeriodicElements = sizeof(elements) / sizeof(PeriodicElement);
166
167static int lineNum = 0;
168
169static const char *
170Itoa(int number)
171{
172    static char buf[200];
173    sprintf(buf, "%d", number);
174    return buf;
175}
176
177static void
178ComputeBonds(Tcl_HashTable *atomTablePtr, Tcl_HashTable *conectTablePtr)
179{
180    PdbAtom **array;
181    int i;
182    Tcl_HashEntry *hPtr;
183    Tcl_HashSearch iter;
184
185#define TOLERANCE 0.45                  /* Fuzz factor for comparing distances
186                                         * (in angstroms) */
187    array = calloc(atomTablePtr->numEntries, sizeof(PdbAtom *));
188    if (array == NULL) {
189        return;
190    }
191    for (i = 0, hPtr = Tcl_FirstHashEntry(atomTablePtr, &iter); hPtr != NULL;
192         hPtr = Tcl_NextHashEntry(&iter), i++) {
193        PdbAtom *atomPtr = Tcl_GetHashValue(hPtr);
194        array[i] = atomPtr;
195    }
196    for (i = 0; i < atomTablePtr->numEntries; i++) {
197        PdbAtom *atom1Ptr;
198        double r1;
199        int j;
200
201        atom1Ptr = array[i];
202        r1 = elements[atom1Ptr->number].radius;
203        for (j = i+1; j < atomTablePtr->numEntries; j++) {
204            PdbAtom *atom2Ptr;
205            double ds2, cut;
206            double r2;
207
208            atom2Ptr  = array[j];
209            if ((atom2Ptr->number == 1) && (atom1Ptr->number == 1)) {
210                continue;
211            }
212            r2 = elements[atom2Ptr->number].radius;
213            cut = (r1 + r2 + TOLERANCE);
214            ds2 = (((atom1Ptr->x - atom2Ptr->x) * (atom1Ptr->x - atom2Ptr->x)) +
215                   ((atom1Ptr->y - atom2Ptr->y) * (atom1Ptr->y - atom2Ptr->y)) +
216                   ((atom1Ptr->z - atom2Ptr->z) * (atom1Ptr->z - atom2Ptr->z)));
217
218            // perform distance test, but ignore pairs between atoms
219            // with nearly identical coords
220            if (ds2 < 0.16) {
221                continue;
222            }
223            if (ds2 < (cut*cut)) {
224                ConnectKey key;
225                int isNew;
226
227                if (atom1Ptr->ordinal > atom2Ptr->ordinal) {
228                    key.from = atom2Ptr->ordinal;
229                    key.to = atom1Ptr->ordinal;
230                } else {
231                    key.from = atom1Ptr->ordinal;
232                    key.to = atom2Ptr->ordinal;
233                }
234                Tcl_CreateHashEntry(conectTablePtr, (char *)&key, &isNew);
235                if (isNew) {
236                    atom1Ptr->numConnections++;
237                    atom2Ptr->numConnections++;
238                }
239            }
240        }
241    }
242    free(array);
243}
244
245int
246GetAtomicNumber(Tcl_Interp *interp, PdbAtom *atomPtr, const char *atomName,
247                const char *symbolName, const char *residueName)
248{
249    char name[5], *namePtr;
250    int elemIndex, symbolIndex;
251    int count;
252    const char *p;
253
254    symbolIndex = elemIndex = -1;
255
256    /*
257     * The atomic symbol may be found the atom name in various locations
258     *  "C   "
259     *  " C  "
260     *  "  C "
261     "  "   C"
262     "  "C 12"
263     "  " C12"
264     */
265    count = 0;
266    for (p = atomName; *p != '\0'; p++) {
267        if (isspace(*p) || isdigit(*p)) {
268            continue;
269        }
270        name[count] = *p;
271        count++;
272    }
273    name[count] = '\0';
274    if (count > 0) {
275        int i;
276
277        namePtr = name;
278        if (name[0] == ' ') {
279            namePtr++;
280        }
281        if (name[1] == ' ') {
282            name[1] = '\0';
283        }
284        for (i = 0; i < numPeriodicElements; i++) {
285            if (strcasecmp(namePtr, elements[i].symbol) == 0) {
286                elemIndex = i;
287                break;
288            }
289        }
290        if ((elemIndex == -1) && (count > 1)) {
291            name[1] = '\0';
292            for (i = 0; i < numPeriodicElements; i++) {
293                if (strcasecmp(namePtr, elements[i].symbol) == 0) {
294                    elemIndex = i;
295                    break;
296                }
297            }
298        }
299    }
300    name[0] = symbolName[0];
301    name[1] = symbolName[1];
302    name[2] = '\0';
303    if (isdigit(name[1])) {
304        sscanf(name, "%d", &atomPtr->number);
305        return TCL_OK;
306    }
307    if ((name[0] != ' ') || (name[1] != ' ')) {
308        int i;
309
310        namePtr = name;
311        if (name[0] == ' ') {
312            namePtr++;
313        }
314        if (name[1] == ' ') {
315            name[1] = '\0';
316        }
317        for (i = 0; i < numPeriodicElements; i++) {
318            if (strcasecmp(namePtr, elements[i].symbol) == 0) {
319                symbolIndex = i;
320                break;
321            }
322        }
323    }
324    if (symbolIndex > 0) {
325        if (elemIndex < 0) {
326            atomPtr->number = symbolIndex;
327            return TCL_OK;
328        }
329        if (symbolIndex != elemIndex) {
330            fprintf(stderr, "atomName %s and symbolName %s don't match\n",
331                    atomName, symbolName);
332            atomPtr->number = symbolIndex;
333            return TCL_OK;
334        }
335        atomPtr->number = elemIndex;
336        return TCL_OK;
337    } else if (elemIndex > 0) {
338        atomPtr->number = elemIndex;
339        return TCL_OK;
340    }
341    Tcl_AppendResult(interp, "line ", Itoa(lineNum), ": bad atom \"",
342                     atomName, "\" or element \"",
343                     symbolName, "\"", (char *)NULL);
344    return TCL_ERROR;
345}
346
347static PdbAtom *
348NewAtom(Tcl_Interp *interp, int ordinal, const char *line, int lineLength)
349{
350    PdbAtom *atomPtr;
351    char buf[200];
352    char atomName[6];                   /* Atom name.  */
353    char symbolName[3];                 /* Symbol name.  */
354    char residueName[4];                /* Residue name. */
355    double value;
356
357    atomPtr = calloc(1, sizeof(PdbAtom));
358    if (atomPtr == NULL) {
359        return NULL;
360    }
361    atomPtr->ordinal = ordinal;
362    strncpy(atomName, line + 12, 4);
363    atomName[4] = '\0';
364    strncpy(residueName, line + 17, 3);
365    residueName[3] = '\0';
366
367    strncpy(buf, line + 30, 8);
368    buf[8] = '\0';
369    if (Tcl_GetDouble(interp, buf, &value) != TCL_OK) {
370        Tcl_AppendResult(interp, "bad x-coordinate \"", buf,
371                         "\"", (char *)NULL);
372        return NULL;
373    }
374    atomPtr->x = value;
375    strncpy(buf, line + 38, 8);
376    buf[8] = '\0';
377    if (Tcl_GetDouble(interp, buf, &value) != TCL_OK) {
378        Tcl_AppendResult(interp, "bad y-coordinate \"", buf,
379                         "\"", (char *)NULL);
380        return NULL;
381    }
382    atomPtr->y = value;
383    strncpy(buf, line + 46, 8);
384    buf[8] = '\0';
385    if (Tcl_GetDouble(interp, buf, &value) != TCL_OK) {
386        Tcl_AppendResult(interp, "bad z-coordinate \"", buf,
387                         "\"", (char *)NULL);
388        return NULL;
389    }
390    atomPtr->z = value;
391    symbolName[2] = symbolName[1]  = symbolName[0] = '\0';
392    if (lineLength >= 77) {
393        symbolName[0] = line[76];
394        symbolName[1] = line[77];
395        symbolName[2] = '\0';
396    }
397    if (GetAtomicNumber(interp, atomPtr, atomName, symbolName, residueName)
398        != TCL_OK) {
399        return NULL;
400    }
401    return atomPtr;
402}
403
404static void
405FreeAtoms(Tcl_HashTable *tablePtr)
406{
407    Tcl_HashEntry *hPtr;
408    Tcl_HashSearch iter;
409
410    for (hPtr = Tcl_FirstHashEntry(tablePtr, &iter); hPtr != NULL;
411         hPtr = Tcl_NextHashEntry(&iter)) {
412        PdbAtom *atomPtr;
413       
414        atomPtr = Tcl_GetHashValue(hPtr);
415        free(atomPtr);
416    }
417    Tcl_DeleteHashTable(tablePtr);
418}
419
420
421static INLINE const char *
422SkipSpaces(const char *string)
423{
424    while (isspace(*string)) {
425        string++;
426    }
427    return string;
428}
429
430static int
431IsSpaces(const char *string)
432{
433    const char *p;
434    for (p = string; *p != '\0'; p++) {
435        if (!isspace(*p)) {
436            return 0;
437        }
438    }
439    return 1;
440}
441
442static INLINE const char *
443GetLine(const char **stringPtr, const char *endPtr, int *lengthPtr)
444{
445    const char *line;
446    const char *p;
447
448    line = SkipSpaces(*stringPtr);
449    for (p = line; p < endPtr; p++) {
450        if (*p == '\n') {
451            *stringPtr = p + 1;
452            *lengthPtr = p - line;
453            return line;
454        }
455    }
456    *lengthPtr = (p - line) - 1;
457    *stringPtr = p;
458    return line;
459}
460
461
462static int
463SerialToAtom(Tcl_Interp *interp, Tcl_HashTable *tablePtr, const char *string,
464             PdbAtom **atomPtrPtr)
465{
466    int serial;
467    long lserial;
468    Tcl_HashEntry *hPtr;
469
470    string = SkipSpaces(string);
471    if (Tcl_GetInt(interp, string, &serial) != TCL_OK) {
472        Tcl_AppendResult(interp, ": line ", Itoa(lineNum),
473                ": invalid serial number \"", string,
474                "\" in CONECT record", (char *)NULL);
475        return TCL_ERROR;
476    }
477    lserial = (long)serial;
478    hPtr = Tcl_FindHashEntry(tablePtr, (char *)lserial);
479    if (hPtr == NULL) {
480        Tcl_AppendResult(interp, "serial number \"", string,
481                         "\" not found in table", (char *)NULL);
482        return TCL_ERROR;
483    }
484    *atomPtrPtr = Tcl_GetHashValue(hPtr);
485    return TCL_OK;
486}
487
488/*
489 *  PdbToVtk string ?-bonds none|both|auto|conect?
490 */
491static int
492PdbToVtkCmd(ClientData clientData, Tcl_Interp *interp, int objc,
493           Tcl_Obj *const *objv)
494{
495    Tcl_Obj *objPtr, *pointsObjPtr, *atomsObjPtr;
496    const char *p, *pend;
497    const char *string;
498    char mesg[2000];
499    int length, nextOrdinal;
500    Tcl_HashTable atomTable, conectTable;
501    Tcl_HashEntry *hPtr;
502    Tcl_HashSearch iter;
503    int isNew;
504    int bondFlags;
505
506    bondFlags = BOND_BOTH;
507    lineNum = nextOrdinal = 0;
508    if ((objc != 2) && (objc != 4)) {
509        Tcl_AppendResult(interp, "wrong # arguments: should be \"",
510                Tcl_GetString(objv[0]),
511                " string ?-bonds none|both|auto|conect\"", (char *)NULL);
512        return TCL_ERROR;
513    }
514    if (objc == 4) {
515        const char *string;
516        char c;
517
518        string = Tcl_GetString(objv[2]);
519        if ((string[0] != '-') || (strcmp(string, "-bonds") != 0)) {
520            Tcl_AppendResult(interp, "bad switch \"", string,
521                             "\": should be \"-bonds\"", (char *)NULL);
522            return TCL_ERROR;
523        }
524        string = Tcl_GetString(objv[3]);
525        c = string[0];
526        if ((c == 'n') && (strcmp(string, "none") == 0)) {
527            bondFlags = BOND_NONE;
528        } else if ((c == 'b') && (strcmp(string, "both") == 0)) {
529            bondFlags = BOND_BOTH;
530        } else if ((c == 'a') && (strcmp(string, "auto") == 0)) {
531            bondFlags = BOND_COMPUTE;
532        } else if ((c == 'c') && (strcmp(string, "conect") == 0)) {
533            bondFlags = BOND_CONECT;
534        } else {
535            Tcl_AppendResult(interp, "bad value \"", string,
536                "\" for -bonds switch: should be none, both, manual, or auto",
537                (char *)NULL);
538            return TCL_ERROR;
539        }
540    }
541    Tcl_InitHashTable(&atomTable, TCL_ONE_WORD_KEYS);
542    Tcl_InitHashTable(&conectTable, sizeof(ConnectKey) / sizeof(int));
543    string = Tcl_GetStringFromObj(objv[1], &length);
544    pointsObjPtr = Tcl_NewStringObj("", -1);
545    atomsObjPtr = Tcl_NewStringObj("", -1);
546    Tcl_IncrRefCount(pointsObjPtr);
547    Tcl_IncrRefCount(atomsObjPtr);
548    objPtr = NULL;
549    for (p = string, pend = p + length; p < pend; /*empty*/) {
550        const char *line;
551        char c;
552        int lineLength;
553
554        lineLength = 0;                 /* Suppress compiler warning. */
555        lineNum++;
556        line = GetLine(&p, pend, &lineLength);
557        if (line >= pend) {
558            break;                      /* EOF */
559        }
560        if ((line[0] == '#') || (line[0] == '\n')) {
561            continue;                   /* Skip blank or comment lines. */
562        }
563        c = line[0];
564        if (((c == 'A') && (strncmp(line, "ATOM  ", 6) == 0)) ||
565            ((c == 'H') && (strncmp(line, "HETATM", 6) == 0))) {
566            char buf[200];
567            int serialNum;
568            long lserial;
569            PdbAtom *atomPtr;
570
571            if (lineLength < 47) {
572                char *msg;
573
574                msg = (char *)line;
575                msg[lineLength] = '\0';
576                Tcl_AppendResult(interp, "line ", Itoa(lineNum),
577                        " short ATOM line \"", line, "\"", (char *)NULL);
578                return TCL_ERROR;
579            }           
580            strncpy(buf, line + 6, 5);
581            buf[5] = '\0';
582            if (Tcl_GetInt(interp, buf, &serialNum) != TCL_OK) {
583                Tcl_AppendResult(interp, "line ", Itoa(lineNum),
584                        ": bad serial number \"", buf, "\"", (char *)NULL);
585                goto error;
586            }
587            lserial = (long)serialNum;
588            hPtr = Tcl_CreateHashEntry(&atomTable, (char *)lserial, &isNew);
589            if (!isNew) {
590                Tcl_AppendResult(interp, "line ", Itoa(lineNum),
591                        ": serial number \"", buf, "\" found more than once",
592                        (char *)NULL);
593                goto error;
594            }
595            atomPtr = NewAtom(interp, nextOrdinal, line, lineLength);
596            if (atomPtr == NULL) {
597                goto error;
598            }
599            Tcl_SetHashValue(hPtr, atomPtr);
600            Tcl_ListObjAppendElement(interp, pointsObjPtr,
601                                     Tcl_NewDoubleObj(atomPtr->x));
602            Tcl_ListObjAppendElement(interp, pointsObjPtr,
603                                     Tcl_NewDoubleObj(atomPtr->y));
604            Tcl_ListObjAppendElement(interp, pointsObjPtr,
605                                     Tcl_NewDoubleObj(atomPtr->z));
606            Tcl_ListObjAppendElement(interp, atomsObjPtr,
607                                     Tcl_NewIntObj(atomPtr->number));
608            nextOrdinal++;
609        } else if ((c == 'C') && (strncmp(line, "CONECT", 6) == 0)) {
610            PdbAtom *atom1Ptr;
611            int i, n;
612            char buf[200];
613
614            if ((bondFlags & BOND_CONECT) == 0) {
615                continue;
616            }
617            strncpy(buf, line + 6, 5);
618            buf[5] = '\0';
619            if (lineLength < 11) {
620                char *msg;
621
622                msg = (char *)line;
623                msg[lineLength] = '\0';
624                Tcl_AppendResult(interp, "line ", Itoa(lineNum),
625                        ": bad CONECT record \"", msg, "\"",
626                        (char *)NULL);
627                goto error;
628            }
629            if (SerialToAtom(interp, &atomTable, buf, &atom1Ptr) != TCL_OK) {
630                goto error;
631            }
632            /* Allow basic maximum of 4 connections. */
633            for (n = 11, i = 0; i < 4; i++, n += 5) {
634                ConnectKey key;
635                PdbAtom *atom2Ptr;
636
637                if (n >= lineLength) {
638                    break;             
639                }
640                strncpy(buf, line + n, 5);
641                buf[5] = '\0';
642                if (IsSpaces(buf)) {
643                    break;              /* No more entries */
644                }
645                if (SerialToAtom(interp, &atomTable, buf, &atom2Ptr) !=TCL_OK) {
646                    goto error;
647                }
648                if (atom1Ptr->ordinal > atom2Ptr->ordinal) {
649                    key.from = atom2Ptr->ordinal;
650                    key.to = atom1Ptr->ordinal;
651                } else {
652                    key.from = atom1Ptr->ordinal;
653                    key.to = atom2Ptr->ordinal;
654                }
655                Tcl_CreateHashEntry(&conectTable, (char *)&key, &isNew);
656                if (isNew) {
657                    atom1Ptr->numConnections++;
658                    atom2Ptr->numConnections++;
659                }
660            }
661        }
662    }
663    if (bondFlags & BOND_COMPUTE) {
664        ComputeBonds(&atomTable, &conectTable);
665    }
666    objPtr = Tcl_NewStringObj("# vtk DataFile Version 2.0\n", -1);
667    Tcl_AppendToObj(objPtr, "Converted from PDB format\n", -1);
668    Tcl_AppendToObj(objPtr, "ASCII\n", -1);
669    Tcl_AppendToObj(objPtr, "DATASET POLYDATA\n", -1);
670    sprintf(mesg, "POINTS %d float\n", atomTable.numEntries);
671    Tcl_AppendToObj(objPtr, mesg, -1);
672    Tcl_AppendObjToObj(objPtr, pointsObjPtr);
673    sprintf(mesg, "\n");
674    Tcl_AppendToObj(objPtr, mesg, -1);
675
676    if (conectTable.numEntries > 0) {
677        sprintf(mesg, "LINES %d %d\n", conectTable.numEntries,
678                conectTable.numEntries * 3);
679        Tcl_AppendToObj(objPtr, mesg, -1);
680        for (hPtr = Tcl_FirstHashEntry(&conectTable, &iter); hPtr != NULL;
681             hPtr = Tcl_NextHashEntry(&iter)) {
682            ConnectKey *keyPtr;
683
684            keyPtr = (ConnectKey *)Tcl_GetHashKey(&conectTable, hPtr);
685            sprintf(mesg, "2 %d %d\n", keyPtr->from, keyPtr->to);
686            Tcl_AppendToObj(objPtr, mesg, -1);
687        }
688        sprintf(mesg, "\n");
689        Tcl_AppendToObj(objPtr, mesg, -1);
690    }
691#if 0
692    for (hPtr = Tcl_FirstHashEntry(&atomTable, &iter); hPtr != NULL;
693         hPtr = Tcl_NextHashEntry(&iter)) {
694        PdbAtom *atomPtr = Tcl_GetHashValue(hPtr);
695        fprintf(stderr, "%d %s %d connections\n", atomPtr->ordinal,
696                elements[atomPtr->number].symbol, atomPtr->numConnections);
697    }
698#endif
699    sprintf(mesg, "POINT_DATA %d\n", atomTable.numEntries);
700    Tcl_AppendToObj(objPtr, mesg, -1);
701    sprintf(mesg, "SCALARS element int\n");
702    Tcl_AppendToObj(objPtr, mesg, -1);
703    sprintf(mesg, "LOOKUP_TABLE default\n");
704    Tcl_AppendToObj(objPtr, mesg, -1);
705    Tcl_AppendObjToObj(objPtr, atomsObjPtr);
706    FreeAtoms(&atomTable);
707    Tcl_DeleteHashTable(&conectTable);
708    Tcl_DecrRefCount(pointsObjPtr);
709    Tcl_DecrRefCount(atomsObjPtr);
710    if (objPtr != NULL) {
711        Tcl_SetObjResult(interp, objPtr);
712    }
713    return TCL_OK;
714 error:
715    Tcl_DeleteHashTable(&conectTable);
716    FreeAtoms(&atomTable);
717    Tcl_DecrRefCount(pointsObjPtr);
718    Tcl_DecrRefCount(atomsObjPtr);
719    return TCL_ERROR;
720}
721
722/*
723 * ------------------------------------------------------------------------
724 *  RpPdbToVtk_Init --
725 *
726 *  Invoked when the Rappture GUI library is being initialized
727 *  to install the "PdbToVtk" command into the interpreter.
728 *
729 *  Returns TCL_OK if successful, or TCL_ERROR (along with an error
730 *  message in the interp) if anything goes wrong.
731 * ------------------------------------------------------------------------
732 */
733int
734RpPdbToVtk_Init(Tcl_Interp *interp)
735{
736    /* install the widget command */
737    Tcl_CreateObjCommand(interp, "Rappture::PdbToVtk", PdbToVtkCmd,
738        NULL, NULL);
739    return TCL_OK;
740}
Note: See TracBrowser for help on using the repository browser.