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

Last change on this file since 3848 was 3848, checked in by gah, 11 years ago

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