source: branches/blt4/gui/src/RpPdbToVtk.c @ 3954

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

sync with trunk
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;
194
195        atomPtr = Tcl_GetHashValue(hPtr);
196        array[i] = atomPtr;
197    }
198    for (i = 0; i < atomTablePtr->numEntries; i++) {
199        PdbAtom *atom1Ptr;
200        double r1;
201        int j;
202
203        atom1Ptr = array[i];
204        r1 = elements[atom1Ptr->number].radius;
205        for (j = i+1; j < atomTablePtr->numEntries; j++) {
206            PdbAtom *atom2Ptr;
207            double ds2, cut;
208            double r2;
209
210            atom2Ptr  = array[j];
211            if ((atom2Ptr->number == 1) && (atom1Ptr->number == 1)) {
212                continue;
213            }
214            r2 = elements[atom2Ptr->number].radius;
215            cut = (r1 + r2 + TOLERANCE);
216            ds2 = (((atom1Ptr->x - atom2Ptr->x) * (atom1Ptr->x - atom2Ptr->x)) +
217                   ((atom1Ptr->y - atom2Ptr->y) * (atom1Ptr->y - atom2Ptr->y)) +
218                   ((atom1Ptr->z - atom2Ptr->z) * (atom1Ptr->z - atom2Ptr->z)));
219
220            // perform distance test, but ignore pairs between atoms
221            // with nearly identical coords
222            if (ds2 < 0.16) {
223                continue;
224            }
225            if (ds2 < (cut*cut)) {
226                ConnectKey key;
227                int isNew;
228
229                if (atom1Ptr->ordinal > atom2Ptr->ordinal) {
230                    key.from = atom2Ptr->ordinal;
231                    key.to = atom1Ptr->ordinal;
232                } else {
233                    key.from = atom1Ptr->ordinal;
234                    key.to = atom2Ptr->ordinal;
235                }
236                Tcl_CreateHashEntry(conectTablePtr, (char *)&key, &isNew);
237                if (isNew) {
238                    atom1Ptr->numConnections++;
239                    atom2Ptr->numConnections++;
240                }
241            }
242        }
243    }
244    free(array);
245}
246
247int
248GetAtomicNumber(Tcl_Interp *interp, PdbAtom *atomPtr, const char *atomName,
249                const char *symbolName, const char *residueName)
250{
251    char name[5], *namePtr;
252    int elemIndex, symbolIndex;
253    int count;
254    const char *p;
255
256    symbolIndex = elemIndex = -1;
257
258    /*
259     * The atomic symbol may be found the atom name in various locations
260     *  "C   "
261     *  " C  "
262     *  "  C "
263     "  "   C"
264     "  "C 12"
265     "  " C12"
266     */
267    count = 0;
268    for (p = atomName; *p != '\0'; p++) {
269        if (isspace(*p) || isdigit(*p)) {
270            continue;
271        }
272        name[count] = *p;
273        count++;
274    }
275    name[count] = '\0';
276    if (count > 0) {
277        int i;
278
279        namePtr = name;
280        if (name[0] == ' ') {
281            namePtr++;
282        }
283        if (name[1] == ' ') {
284            name[1] = '\0';
285        }
286        for (i = 0; i < numPeriodicElements; i++) {
287            if (strcasecmp(namePtr, elements[i].symbol) == 0) {
288                elemIndex = i;
289                break;
290            }
291        }
292        if ((elemIndex == -1) && (count > 1)) {
293            name[1] = '\0';
294            for (i = 0; i < numPeriodicElements; i++) {
295                if (strcasecmp(namePtr, elements[i].symbol) == 0) {
296                    elemIndex = i;
297                    break;
298                }
299            }
300        }
301    }
302    name[0] = symbolName[0];
303    name[1] = symbolName[1];
304    name[2] = '\0';
305    if (isdigit(name[1])) {
306        sscanf(name, "%d", &atomPtr->number);
307        return TCL_OK;
308    }
309    if ((name[0] != ' ') || (name[1] != ' ')) {
310        int i;
311
312        namePtr = name;
313        if (name[0] == ' ') {
314            namePtr++;
315        }
316        if (name[1] == ' ') {
317            name[1] = '\0';
318        }
319        for (i = 0; i < numPeriodicElements; i++) {
320            if (strcasecmp(namePtr, elements[i].symbol) == 0) {
321                symbolIndex = i;
322                break;
323            }
324        }
325    }
326    if (symbolIndex > 0) {
327        if (elemIndex < 0) {
328            atomPtr->number = symbolIndex;
329            return TCL_OK;
330        }
331        if (symbolIndex != elemIndex) {
332            fprintf(stderr, "atomName %s and symbolName %s don't match\n",
333                    atomName, symbolName);
334            atomPtr->number = symbolIndex;
335            return TCL_OK;
336        }
337        atomPtr->number = elemIndex;
338        return TCL_OK;
339    } else if (elemIndex > 0) {
340        atomPtr->number = elemIndex;
341        return TCL_OK;
342    }
343    Tcl_AppendResult(interp, "line ", Itoa(lineNum), ": bad atom \"",
344                     atomName, "\" or element \"",
345                     symbolName, "\"", (char *)NULL);
346    return TCL_ERROR;
347}
348
349static PdbAtom *
350NewAtom(Tcl_Interp *interp, int ordinal, const char *line, int lineLength)
351{
352    PdbAtom *atomPtr;
353    char buf[200];
354    char atomName[6];                   /* Atom name.  */
355    char symbolName[3];                 /* Symbol name.  */
356    char residueName[4];                /* Residue name. */
357    double value;
358
359    atomPtr = calloc(1, sizeof(PdbAtom));
360    if (atomPtr == NULL) {
361        return NULL;
362    }
363    atomPtr->ordinal = ordinal;
364    strncpy(atomName, line + 12, 4);
365    atomName[4] = '\0';
366    strncpy(residueName, line + 17, 3);
367    residueName[3] = '\0';
368
369    strncpy(buf, line + 30, 8);
370    buf[8] = '\0';
371    if (Tcl_GetDouble(interp, buf, &value) != TCL_OK) {
372        Tcl_AppendResult(interp, "bad x-coordinate \"", buf,
373                         "\"", (char *)NULL);
374        return NULL;
375    }
376    atomPtr->x = value;
377    strncpy(buf, line + 38, 8);
378    buf[8] = '\0';
379    if (Tcl_GetDouble(interp, buf, &value) != TCL_OK) {
380        Tcl_AppendResult(interp, "bad y-coordinate \"", buf,
381                         "\"", (char *)NULL);
382        return NULL;
383    }
384    atomPtr->y = value;
385    strncpy(buf, line + 46, 8);
386    buf[8] = '\0';
387    if (Tcl_GetDouble(interp, buf, &value) != TCL_OK) {
388        Tcl_AppendResult(interp, "bad z-coordinate \"", buf,
389                         "\"", (char *)NULL);
390        return NULL;
391    }
392    atomPtr->z = value;
393    symbolName[2] = symbolName[1]  = symbolName[0] = '\0';
394    if (lineLength >= 77) {
395        symbolName[0] = line[76];
396        symbolName[1] = line[77];
397        symbolName[2] = '\0';
398    }
399    if (GetAtomicNumber(interp, atomPtr, atomName, symbolName, residueName)
400        != TCL_OK) {
401        return NULL;
402    }
403    return atomPtr;
404}
405
406static void
407FreeAtoms(Tcl_HashTable *tablePtr)
408{
409    Tcl_HashEntry *hPtr;
410    Tcl_HashSearch iter;
411
412    for (hPtr = Tcl_FirstHashEntry(tablePtr, &iter); hPtr != NULL;
413         hPtr = Tcl_NextHashEntry(&iter)) {
414        PdbAtom *atomPtr;
415       
416        atomPtr = Tcl_GetHashValue(hPtr);
417        free(atomPtr);
418    }
419    Tcl_DeleteHashTable(tablePtr);
420}
421
422
423static INLINE const char *
424SkipSpaces(const char *string)
425{
426    while (isspace(*string)) {
427        string++;
428    }
429    return string;
430}
431
432static int
433IsSpaces(const char *string)
434{
435    const char *p;
436    for (p = string; *p != '\0'; p++) {
437        if (!isspace(*p)) {
438            return 0;
439        }
440    }
441    return 1;
442}
443
444static INLINE const char *
445GetLine(const char **stringPtr, const char *endPtr, int *lengthPtr)
446{
447    const char *line;
448    const char *p;
449
450    line = SkipSpaces(*stringPtr);
451    for (p = line; p < endPtr; p++) {
452        if (*p == '\n') {
453            *stringPtr = p + 1;
454            *lengthPtr = p - line;
455            return line;
456        }
457    }
458    *lengthPtr = (p - line) - 1;
459    *stringPtr = p;
460    return line;
461}
462
463
464static int
465SerialToAtom(Tcl_Interp *interp, Tcl_HashTable *tablePtr, const char *string,
466                PdbAtom **atomPtrPtr)
467{
468    int serial;
469    long lserial;
470    Tcl_HashEntry *hPtr;
471
472    string = SkipSpaces(string);
473    if (Tcl_GetInt(interp, string, &serial) != TCL_OK) {
474        Tcl_AppendResult(interp, ": line ", Itoa(lineNum),
475                ": invalid serial number \"", string,
476                "\" in CONECT record", (char *)NULL);
477        return TCL_ERROR;
478    }
479    lserial = (long)serial;
480    hPtr = Tcl_FindHashEntry(tablePtr, (char *)lserial);
481    if (hPtr == NULL) {
482        Tcl_AppendResult(interp, "serial number \"", string,
483                         "\" not found in table", (char *)NULL);
484        return TCL_ERROR;
485    }
486    *atomPtrPtr = Tcl_GetHashValue(hPtr);
487    return TCL_OK;
488}
489
490/*
491 *  PdbToVtk string ?-bonds none|both|auto|conect?
492 */
493static int
494PdbToVtkCmd(ClientData clientData, Tcl_Interp *interp, int objc,
495           Tcl_Obj *const *objv)
496{
497    Tcl_Obj *objPtr, *pointsObjPtr, *atomsObjPtr;
498    const char *p, *pend;
499    const char *string;
500    char mesg[2000];
501    int length, nextOrdinal;
502    Tcl_HashTable atomTable, conectTable;
503    Tcl_HashEntry *hPtr;
504    Tcl_HashSearch iter;
505    int isNew;
506    int bondFlags;
507
508    bondFlags = BOND_BOTH;
509    lineNum = nextOrdinal = 0;
510    if ((objc != 2) && (objc != 4)) {
511        Tcl_AppendResult(interp, "wrong # arguments: should be \"",
512                Tcl_GetString(objv[0]),
513                " string ?-bonds none|both|auto|conect\"", (char *)NULL);
514        return TCL_ERROR;
515    }
516    if (objc == 4) {
517        const char *string;
518        char c;
519
520        string = Tcl_GetString(objv[2]);
521        if ((string[0] != '-') || (strcmp(string, "-bonds") != 0)) {
522            Tcl_AppendResult(interp, "bad switch \"", string,
523                             "\": should be \"-bonds\"", (char *)NULL);
524            return TCL_ERROR;
525        }
526        string = Tcl_GetString(objv[3]);
527        c = string[0];
528        if ((c == 'n') && (strcmp(string, "none") == 0)) {
529            bondFlags = BOND_NONE;
530        } else if ((c == 'b') && (strcmp(string, "both") == 0)) {
531            bondFlags = BOND_BOTH;
532        } else if ((c == 'a') && (strcmp(string, "auto") == 0)) {
533            bondFlags = BOND_COMPUTE;
534        } else if ((c == 'c') && (strcmp(string, "conect") == 0)) {
535            bondFlags = BOND_CONECT;
536        } else {
537            Tcl_AppendResult(interp, "bad value \"", string,
538                "\" for -bonds switch: should be none, both, manual, or auto",
539                (char *)NULL);
540            return TCL_ERROR;
541        }
542    }
543    Tcl_InitHashTable(&atomTable, TCL_ONE_WORD_KEYS);
544    Tcl_InitHashTable(&conectTable, sizeof(ConnectKey) / sizeof(int));
545    string = Tcl_GetStringFromObj(objv[1], &length);
546    pointsObjPtr = Tcl_NewStringObj("", -1);
547    atomsObjPtr = Tcl_NewStringObj("", -1);
548    Tcl_IncrRefCount(pointsObjPtr);
549    Tcl_IncrRefCount(atomsObjPtr);
550    objPtr = NULL;
551    for (p = string, pend = p + length; p < pend; /*empty*/) {
552        const char *line;
553        char c;
554        int lineLength;
555
556        lineLength = 0;                 /* Suppress compiler warning. */
557        lineNum++;
558        line = GetLine(&p, pend, &lineLength);
559        if (line >= pend) {
560            break;                      /* EOF */
561        }
562        if ((line[0] == '#') || (line[0] == '\n')) {
563            continue;                   /* Skip blank or comment lines. */
564        }
565        c = line[0];
566        if (((c == 'A') && (strncmp(line, "ATOM  ", 6) == 0)) ||
567            ((c == 'H') && (strncmp(line, "HETATM", 6) == 0))) {
568            char buf[200];
569            int serialNum;
570            long lserial;
571            PdbAtom *atomPtr;
572
573            if (lineLength < 47) {
574                char *msg;
575
576                msg = (char *)line;
577                msg[lineLength] = '\0';
578                Tcl_AppendResult(interp, "line ", Itoa(lineNum),
579                        " short ATOM line \"", line, "\"", (char *)NULL);
580                return TCL_ERROR;
581            }           
582            strncpy(buf, line + 6, 5);
583            buf[5] = '\0';
584            if (Tcl_GetInt(interp, buf, &serialNum) != TCL_OK) {
585                Tcl_AppendResult(interp, "line ", Itoa(lineNum),
586                        ": bad serial number \"", buf, "\"", (char *)NULL);
587                goto error;
588            }
589            lserial = (long)serialNum;
590            hPtr = Tcl_CreateHashEntry(&atomTable, (char *)lserial, &isNew);
591            if (!isNew) {
592                Tcl_AppendResult(interp, "line ", Itoa(lineNum),
593                        ": serial number \"", buf, "\" found more than once",
594                        (char *)NULL);
595                goto error;
596            }
597            atomPtr = NewAtom(interp, nextOrdinal, line, lineLength);
598            if (atomPtr == NULL) {
599                goto error;
600            }
601            Tcl_SetHashValue(hPtr, atomPtr);
602            Tcl_ListObjAppendElement(interp, pointsObjPtr,
603                                     Tcl_NewDoubleObj(atomPtr->x));
604            Tcl_ListObjAppendElement(interp, pointsObjPtr,
605                                     Tcl_NewDoubleObj(atomPtr->y));
606            Tcl_ListObjAppendElement(interp, pointsObjPtr,
607                                     Tcl_NewDoubleObj(atomPtr->z));
608            Tcl_ListObjAppendElement(interp, atomsObjPtr,
609                                     Tcl_NewIntObj(atomPtr->number));
610            nextOrdinal++;
611        } else if ((c == 'C') && (strncmp(line, "CONECT", 6) == 0)) {
612            PdbAtom *atom1Ptr;
613            int i, n;
614            char buf[200];
615
616            if ((bondFlags & BOND_CONECT) == 0) {
617                continue;
618            }
619            strncpy(buf, line + 6, 5);
620            buf[5] = '\0';
621            if (lineLength < 11) {
622                char *msg;
623
624                msg = (char *)line;
625                msg[lineLength] = '\0';
626                Tcl_AppendResult(interp, "line ", Itoa(lineNum),
627                        ": bad CONECT record \"", msg, "\"",
628                        (char *)NULL);
629                goto error;
630            }
631            if (SerialToAtom(interp, &atomTable, buf, &atom1Ptr) != TCL_OK) {
632                goto error;
633            }
634            /* Allow basic maximum of 4 connections. */
635            for (n = 11, i = 0; i < 4; i++, n += 5) {
636                ConnectKey key;
637                PdbAtom *atom2Ptr;
638
639                if (n >= lineLength) {
640                    break;             
641                }
642                strncpy(buf, line + n, 5);
643                buf[5] = '\0';
644                if (IsSpaces(buf)) {
645                    break;              /* No more entries */
646                }
647                if (SerialToAtom(interp, &atomTable, buf, &atom2Ptr) !=TCL_OK) {
648                    goto error;
649                }
650                if (atom1Ptr->ordinal > atom2Ptr->ordinal) {
651                    key.from = atom2Ptr->ordinal;
652                    key.to = atom1Ptr->ordinal;
653                } else {
654                    key.from = atom1Ptr->ordinal;
655                    key.to = atom2Ptr->ordinal;
656                }
657                Tcl_CreateHashEntry(&conectTable, (char *)&key, &isNew);
658                if (isNew) {
659                    atom1Ptr->numConnections++;
660                    atom2Ptr->numConnections++;
661                }
662            }
663        }
664    }
665    if (bondFlags & BOND_COMPUTE) {
666        ComputeBonds(&atomTable, &conectTable);
667    }
668    objPtr = Tcl_NewStringObj("# vtk DataFile Version 2.0\n", -1);
669    Tcl_AppendToObj(objPtr, "Converted from PDB format\n", -1);
670    Tcl_AppendToObj(objPtr, "ASCII\n", -1);
671    Tcl_AppendToObj(objPtr, "DATASET POLYDATA\n", -1);
672    sprintf(mesg, "POINTS %d float\n", atomTable.numEntries);
673    Tcl_AppendToObj(objPtr, mesg, -1);
674    Tcl_AppendObjToObj(objPtr, pointsObjPtr);
675    sprintf(mesg, "\n");
676    Tcl_AppendToObj(objPtr, mesg, -1);
677
678    if (conectTable.numEntries > 0) {
679        sprintf(mesg, "LINES %d %d\n", conectTable.numEntries,
680                conectTable.numEntries * 3);
681        Tcl_AppendToObj(objPtr, mesg, -1);
682        for (hPtr = Tcl_FirstHashEntry(&conectTable, &iter); hPtr != NULL;
683             hPtr = Tcl_NextHashEntry(&iter)) {
684            ConnectKey *keyPtr;
685
686            keyPtr = (ConnectKey *)Tcl_GetHashKey(&conectTable, hPtr);
687            sprintf(mesg, "2 %d %d\n", keyPtr->from, keyPtr->to);
688            Tcl_AppendToObj(objPtr, mesg, -1);
689        }
690        sprintf(mesg, "\n");
691        Tcl_AppendToObj(objPtr, mesg, -1);
692    }
693    for (hPtr = Tcl_FirstHashEntry(&atomTable, &iter); hPtr != NULL;
694         hPtr = Tcl_NextHashEntry(&iter)) {
695        PdbAtom *atomPtr;
696
697        atomPtr = Tcl_GetHashValue(hPtr);
698        fprintf(stderr, "%d %s %d connections\n", atomPtr->ordinal,
699                elements[atomPtr->number].symbol, atomPtr->numConnections);
700    }
701    sprintf(mesg, "POINT_DATA %d\n", atomTable.numEntries);
702    Tcl_AppendToObj(objPtr, mesg, -1);
703    sprintf(mesg, "SCALARS element int\n");
704    Tcl_AppendToObj(objPtr, mesg, -1);
705    sprintf(mesg, "LOOKUP_TABLE default\n");
706    Tcl_AppendToObj(objPtr, mesg, -1);
707    Tcl_AppendObjToObj(objPtr, atomsObjPtr);
708    FreeAtoms(&atomTable);
709    Tcl_DeleteHashTable(&conectTable);
710    Tcl_DecrRefCount(pointsObjPtr);
711    Tcl_DecrRefCount(atomsObjPtr);
712    if (objPtr != NULL) {
713        Tcl_SetObjResult(interp, objPtr);
714    }
715    return TCL_OK;
716 error:
717    Tcl_DeleteHashTable(&conectTable);
718    FreeAtoms(&atomTable);
719    Tcl_DecrRefCount(pointsObjPtr);
720    Tcl_DecrRefCount(atomsObjPtr);
721    return TCL_ERROR;
722}
723
724/*
725 * ------------------------------------------------------------------------
726 *  RpPdbToVtk_Init --
727 *
728 *  Invoked when the Rappture GUI library is being initialized
729 *  to install the "PdbToVtk" command into the interpreter.
730 *
731 *  Returns TCL_OK if successful, or TCL_ERROR (along with an error
732 *  message in the interp) if anything goes wrong.
733 * ------------------------------------------------------------------------
734 */
735int
736RpPdbToVtk_Init(Tcl_Interp *interp)
737{
738    /* install the widget command */
739    Tcl_CreateObjCommand(interp, "Rappture::PdbToVtk", PdbToVtkCmd,
740        NULL, NULL);
741    return TCL_OK;
742}
Note: See TracBrowser for help on using the repository browser.