ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-3.0/src/openbabel/chains.cpp
Revision: 2440
Committed: Wed Nov 16 19:42:11 2005 UTC (18 years, 7 months ago) by tim
File size: 63426 byte(s)
Log Message:
adding openbabel

File Contents

# User Rev Content
1 tim 2440 /**********************************************************************
2     chains.cpp - Parse for macromolecule chains and residues.
3    
4     Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
5     (original author, Roger Sayle, version 1.6, March 1998)
6     (modified by Joe Corkery, OpenEye Scientific Software, March 2001)
7     Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
8    
9     This file is part of the Open Babel project.
10     For more information, see <http://openbabel.sourceforge.net/>
11    
12     This program is free software; you can redistribute it and/or modify
13     it under the terms of the GNU General Public License as published by
14     the Free Software Foundation version 2 of the License.
15    
16     This program is distributed in the hope that it will be useful,
17     but WITHOUT ANY WARRANTY; without even the implied warranty of
18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19     GNU General Public License for more details.
20     ***********************************************************************/
21    
22     //////////////////////////////////////////////////////////////////////////////
23     // File Includes
24     //////////////////////////////////////////////////////////////////////////////
25    
26     #include <stdlib.h>
27     #include <string.h>
28     #include <stdio.h>
29     #include <ctype.h>
30     #include <map>
31    
32     #include "mol.hpp"
33     #include "chains.hpp"
34    
35     using namespace std;
36    
37     //////////////////////////////////////////////////////////////////////////////
38     // Preprocessor Definitions
39     //////////////////////////////////////////////////////////////////////////////
40    
41     #define RESIDMIN 3
42     #define RESIDMAX 32
43     #define ATOMMAX 68
44     #define ATOMMINAMINO 4
45     #define ATOMMINNUCLEIC 50
46     #define ELEMMAX 104
47     #define MAXPEPTIDE 11
48     #define MAXNUCLEIC 15
49     #define AMINOMAX 21
50     #define NUCLEOMAX 6
51     #define STACKSIZE 20
52     #define BUFMAX 8192
53     #define MAXCOVAL 2.0
54     #define SLOPFACTOR 0.56
55     #define THRESHOLD 12
56    
57     #define AI_N 0
58     #define AI_CA 1
59     #define AI_C 2
60     #define AI_O 3
61     #define AI_OXT 37
62    
63     #define AI_P 38
64     #define AI_O1P 39
65     #define AI_O2P 40
66     #define AI_O5 41
67     #define AI_C5 42
68     #define AI_C4 43
69     #define AI_O4 44
70     #define AI_C3 45
71     #define AI_O3 46
72     #define AI_C2 47
73     #define AI_O2 48
74     #define AI_C1 49
75    
76     #define BitN 0x0001
77     #define BitNTer 0x0002
78     #define BitNPro 0x0004
79     #define BitNPT 0x0008
80     #define BitCA 0x0010
81     #define BitCAGly 0x0020
82     #define BitC 0x0100
83     #define BitCTer 0x0200
84     #define BitCOXT 0x0400
85     #define BitO 0x1000
86     #define BitOXT 0x2000
87    
88     #define BitNAll 0x000F
89     #define BitCAAll 0x0030
90     #define BitCAll 0x0700
91     #define BitOAll 0x3000
92    
93     #define BitP 0x0001
94     #define BitPTer 0x0002
95     #define BitOP 0x0004
96     #define BitO5 0x0008
97     #define BitO5Ter 0x0010
98     #define BitC5 0x0020
99     #define BitC4 0x0040
100     #define BitO4 0x0080
101     #define BitC3 0x0100
102     #define BitO3 0x0200
103     #define BitO3Ter 0x0400
104     #define BitC2RNA 0x0800
105     #define BitC2DNA 0x1000
106     #define BitO2 0x2000
107     #define BitC1 0x4000
108    
109     #define BitPAll 0x0003
110     #define Bit05All 0x0018
111     #define BitO3All 0x0600
112     #define BitC2All 0x1800
113    
114     #define BC_ASSIGN 0x01
115     #define BC_COUNT 0x02
116     #define BC_ELEM 0x03
117     #define BC_EVAL 0x04
118     #define BC_IDENT 0x05
119     #define BC_LOCAL 0x06
120    
121     #define BF_SINGLE 0x01
122     #define BF_DOUBLE 0x02
123     #define BF_TRIPLE 0x04
124     #define BF_AROMATIC 0x08
125    
126     namespace OpenBabel
127     {
128    
129     OBChainsParser chainsparser;
130    
131     //////////////////////////////////////////////////////////////////////////////
132     // Structure / Type Definitions
133     //////////////////////////////////////////////////////////////////////////////
134    
135     typedef struct
136     {
137     char *name;
138     char *data;
139     }
140     ResidType;
141    
142     typedef struct
143     {
144     int atomid,elem;
145     int bcount;
146     int index;
147     }
148     MonoAtomType;
149    
150     typedef struct
151     {
152     int src,dst;
153     int index;
154     int flag;
155     }
156     MonoBondType;
157    
158     typedef struct
159     {
160     int type;
161     union _ByteCode *next;
162     }
163     MonOpStruct;
164    
165     typedef struct
166     {
167     int type;
168     int value;
169     union _ByteCode *tcond;
170     union _ByteCode *fcond;
171     }
172     BinOpStruct;
173    
174     typedef struct
175     {
176     int type;
177     int resid;
178     int *atomid;
179     int *bflags;
180     }
181     AssignStruct;
182    
183     typedef union _ByteCode
184     {
185     int type;
186     MonOpStruct eval; /* BC_EVAL */
187     BinOpStruct count; /* BC_COUNT */
188     BinOpStruct elem; /* BC_ELEM */
189     BinOpStruct ident; /* BC_IDENT */
190     BinOpStruct local; /* BC_LOCAL */
191     AssignStruct assign; /* BC_ASSIGN */
192     } ByteCode;
193    
194     typedef struct
195     {
196     int atom,bond;
197     int prev;
198     }
199     StackType;
200    
201     //////////////////////////////////////////////////////////////////////////////
202     // Global Variables / Tables
203     //////////////////////////////////////////////////////////////////////////////
204    
205     static char ChainsAtomName[ATOMMAX][4] = {
206     /* 0 */ { ' ', 'N', ' ', ' ' },
207     /* 1 */ { ' ', 'C', 'A', ' ' },
208     /* 2 */ { ' ', 'C', ' ', ' ' },
209     /* 3 */ { ' ', 'O', ' ', ' ' },
210     /* 4 */ { ' ', 'C', 'B', ' ' },
211     /* 5 */ { ' ', 'S', 'G', ' ' },
212     /* 6 */ { ' ', 'O', 'G', ' ' },
213     /* 7 */ { ' ', 'C', 'G', ' ' },
214     /* 8 */ { ' ', 'O', 'G', '1' },
215     /* 9 */ { ' ', 'C', 'G', '1' },
216     /* 10 */ { ' ', 'C', 'G', '2' },
217     /* 11 */ { ' ', 'C', 'D', ' ' },
218     /* 12 */ { ' ', 'O', 'D', ' ' },
219     /* 13 */ { ' ', 'S', 'D', ' ' },
220     /* 14 */ { ' ', 'C', 'D', '1' },
221     /* 15 */ { ' ', 'O', 'D', '1' },
222     /* 16 */ { ' ', 'N', 'D', '1' },
223     /* 17 */ { ' ', 'C', 'D', '2' },
224     /* 18 */ { ' ', 'O', 'D', '2' },
225     /* 19 */ { ' ', 'N', 'D', '2' },
226     /* 20 */ { ' ', 'C', 'E', ' ' },
227     /* 21 */ { ' ', 'N', 'E', ' ' },
228     /* 22 */ { ' ', 'C', 'E', '1' },
229     /* 23 */ { ' ', 'O', 'E', '1' },
230     /* 24 */ { ' ', 'N', 'E', '1' },
231     /* 25 */ { ' ', 'C', 'E', '2' },
232     /* 26 */ { ' ', 'O', 'E', '2' },
233     /* 27 */ { ' ', 'N', 'E', '2' },
234     /* 28 */ { ' ', 'C', 'E', '3' },
235     /* 29 */ { ' ', 'C', 'Z', ' ' },
236     /* 30 */ { ' ', 'N', 'Z', ' ' },
237     /* 31 */ { ' ', 'C', 'Z', '2' },
238     /* 32 */ { ' ', 'C', 'Z', '3' },
239     /* 33 */ { ' ', 'O', 'H', ' ' },
240     /* 34 */ { ' ', 'N', 'H', '1' },
241     /* 35 */ { ' ', 'N', 'H', '2' },
242     /* 36 */ { ' ', 'C', 'H', '2' },
243     /* 37 */ { ' ', 'O', 'X', 'T' },
244    
245     /* 38 */ { ' ', 'P', ' ', ' ' },
246     /* 39 */ { ' ', 'O', '1', 'P' },
247     /* 40 */ { ' ', 'O', '2', 'P' },
248     /* 41 */ { ' ', 'O', '5', '*' },
249     /* 42 */ { ' ', 'C', '5', '*' },
250     /* 43 */ { ' ', 'C', '4', '*' },
251     /* 44 */ { ' ', 'O', '4', '*' },
252     /* 45 */ { ' ', 'C', '3', '*' },
253     /* 46 */ { ' ', 'O', '3', '*' },
254     /* 47 */ { ' ', 'C', '2', '*' },
255     /* 48 */ { ' ', 'O', '2', '*' },
256     /* 49 */ { ' ', 'C', '1', '*' },
257     /* 50 */ { ' ', 'N', '9', ' ' },
258     /* 51 */ { ' ', 'C', '8', ' ' },
259     /* 52 */ { ' ', 'N', '7', ' ' },
260     /* 53 */ { ' ', 'C', '5', ' ' },
261     /* 54 */ { ' ', 'C', '6', ' ' },
262     /* 55 */ { ' ', 'O', '6', ' ' },
263     /* 56 */ { ' ', 'N', '6', ' ' },
264     /* 57 */ { ' ', 'N', '1', ' ' },
265     /* 58 */ { ' ', 'C', '2', ' ' },
266     /* 59 */ { ' ', 'O', '2', ' ' },
267     /* 60 */ { ' ', 'N', '2', ' ' },
268     /* 61 */ { ' ', 'N', '3', ' ' },
269     /* 62 */ { ' ', 'C', '4', ' ' },
270     /* 63 */ { ' ', 'O', '4', ' ' },
271     /* 64 */ { ' ', 'N', '4', ' ' },
272     /* 65 */ { ' ', 'C', '5', ' ' },
273     /* 66 */ { ' ', 'C', '5', 'M' },
274     /* 67 */ { ' ', 'C', '6', ' ' }
275     };
276    
277     static Template Peptide[MAXPEPTIDE] = {
278     /* N */ { 0x0001, 7, 2, 0x0030, 0x0100, 0, 0 },
279     /* NTer */ { 0x0002, 7, 1, 0x0030, 0, 0, 0 },
280     /* NPro */ { 0x0004, 7, 3, 0x0030, 0x0100, -6, 0 },
281     /* NPT */ { 0x0008, 7, 2, 0x0030, -6, 0, 0 },
282     /* CA */ { 0x0010, 6, 3, 0x000F, 0x0700, -6, 0 },
283     /* CAGly */ { 0x0020, 6, 2, 0x0003, 0x0700, 0, 0 },
284     /* C */ { 0x0100, 6, 3, 0x0030, 0x1000, 0x0005, 0 },
285     /* CTer */ { 0x0200, 6, 2, 0x0030, 0x1000, 0, 0 },
286     /* COXT */ { 0x0400, 6, 3, 0x0030, 0x1000, 0x2000, 0 },
287     /* O */ { 0x1000, 8, 1, 0x0700, 0, 0, 0 },
288     /* OXT */ { 0x2000, 8, 1, 0x0400, 0, 0, 0 }
289     };
290    
291     static Template Nucleotide[MAXNUCLEIC] = {
292     /* P */ { 0x0001, 15, 4, 0x0004, 0x0004, 0x0008, 0x0200 },
293     /* PTer */ { 0x0002, 15, 3, 0x0004, 0x0004, 0x0008, 0 },
294     /* OP */ { 0x0004, 8, 1, 0x0003, 0, 0, 0 },
295     /* O5 */ { 0x0008, 8, 2, 0x0020, 0x0003, 0, 0 },
296     /* O5Ter */ { 0x0010, 8, 1, 0x0020, 0, 0, 0 },
297     /* C5 */ { 0x0020, 6, 2, 0x0018, 0x0040, 0, 0 },
298     /* C4 */ { 0x0040, 6, 3, 0x0020, 0x0080, 0x0100, 0 },
299     /* O4 */ { 0x0080, 8, 2, 0x0040, 0x4000, 0, 0 },
300     /* C3 */ { 0x0100, 6, 3, 0x0040, 0x0600, 0x1800, 0 },
301     /* O3 */ { 0x0200, 8, 2, 0x0100, 0x0001, 0, 0 },
302     /* O3Ter */ { 0x0400, 8, 1, 0x0100, 0, 0, 0 },
303     /* C2RNA */ { 0x0800, 6, 3, 0x0100, 0x4000, 0x2000, 0 },
304     /* C2DNA */ { 0x1000, 6, 2, 0x0100, 0x4000, 0, 0 },
305     /* O2 */ { 0x2000, 8, 1, 0x0800, 0, 0, 0 },
306     /* C1 */ { 0x4000, 6, 3, 0x0080, 0x1800, -7, 0 }
307     };
308    
309     static char ChainsResName[RESIDMAX][4] = {
310     /*0*/ "UNK", /*1*/ "HOH", /*2*/ "LIG"
311     };
312    
313     static ResidType AminoAcids[AMINOMAX] = {
314     { "ILE", "1-4(-9-14)-10" },
315     { "VAL", "1-4(-9)-10" },
316    
317     { "ALA", "1-4" },
318     { "ASN", "1-4-7(=15)-19" },
319     { "ASP", "1-4-7(=15)-18" },
320     { "ARG", "1-4-7-11-21-29(=34)-35" },
321     { "CYS", "1-4-5" },
322     { "GLN", "1-4-7-11(=23)-27" },
323     { "GLU", "1-4-7-11(=23)-26" },
324     { "GLY", "1" },
325     { "HIS", "1-4-7^16~22^27^17~7" },
326     { "HYP", "1-4-7(-12)-11-0" }, /* ??? */
327     { "LEU", "1-4-7(-14)-17" },
328     { "LYS", "1-4-7-11-20-30" },
329     { "MET", "1-4-7-13-20" },
330     { "PHE", "1-4-7~14^22~29^25~17^7" },
331     { "PRO", "1-4-7-11-0" },
332     { "SER", "1-4-6" },
333     { "THR", "1-4(-8)-10" },
334     { "TRP", "1-4-7~14^24^25~17(^7)^28~32^36~31^25" },
335     { "TYR", "1-4-7~14^22~29(-33)^25~17^7" }
336     };
337    
338     /* Pyroglutamate (PCA): 1-4-7-11(=" OE ")-0 PDB Example: 1CEL */
339     /* Amino-N-Butyric Acid (ABA): 1-4-7 PDB Example: 1BBO */
340     /* Selenic Acid (SEC): 1-4-"SEG "(-15)-18 PDB Example: 1GP1 */
341    
342     static ResidType Nucleotides[NUCLEOMAX] = {
343     { " A", "49-50-51-52-53-54(-56)-57-58-61-62(-53)-50" },
344     { " C", "49-57-58(-59)-61-62(-64)-65-67-57" },
345     { " G", "49-50-51-52-53-54(-55)-57-58(-60)-61-62(-53)-50" },
346     { " T", "49-57-58(-59)-61-62(-63)-65(-66)-67-57" },
347     { " U", "49-57-58(-59)-61-62(-63)-65-67-57" },
348     { " I", "49-50-51-52-53-54(-55)-57-58-61-62(-53)-50" }
349     };
350    
351     static MonoAtomType MonoAtom[MaxMonoAtom];
352     static MonoBondType MonoBond[MaxMonoBond];
353     static int MonoAtomCount;
354     static int MonoBondCount;
355    
356     static StackType Stack[STACKSIZE];
357     static int StackPtr;
358    
359     static int AtomIndex;
360     static int BondIndex;
361     static bool StrictFlag = false;
362    
363     //////////////////////////////////////////////////////////////////////////////
364     // Static Functions
365     //////////////////////////////////////////////////////////////////////////////
366    
367     static ByteCode *AllocateByteCode(int type)
368     {
369     ByteCode *result;
370    
371     result = (ByteCode*)malloc(sizeof(ByteCode));
372     if( !result )
373     {
374     obErrorLog.ThrowError(__FUNCTION__, "Unable to allocate byte codes for biomolecule residue perception.", obError);
375     // exit(1);
376     }
377     result->type = type;
378     result->eval.next = NULL;
379     result->count.tcond = NULL;
380     result->count.fcond = NULL;
381     result->elem.tcond = NULL;
382     result->elem.fcond = NULL;
383     result->ident.tcond = NULL;
384     result->ident.fcond = NULL;
385     result->local.tcond = NULL;
386     result->local.fcond = NULL;
387     result->assign.atomid = NULL;
388     result->assign.bflags = NULL;
389    
390     return (result);
391     }
392    
393     static void DeleteByteCode(ByteCode *node)
394     {
395     if (node == NULL)
396     return;
397     else
398     {
399     switch (node->type)
400     {
401     case BC_ASSIGN:
402    
403     if (node->assign.atomid != NULL)
404     free(node->assign.atomid);
405     if (node->assign.bflags != NULL)
406     free(node->assign.bflags);
407    
408     break;
409    
410     case BC_COUNT:
411    
412     DeleteByteCode(node->count.tcond);
413     DeleteByteCode(node->count.fcond);
414     break;
415     case BC_ELEM:
416    
417     DeleteByteCode(node->elem.tcond);
418     DeleteByteCode(node->elem.fcond);
419     break;
420    
421     case BC_EVAL:
422    
423     DeleteByteCode(node->eval.next);
424     break;
425    
426     case BC_IDENT:
427    
428     DeleteByteCode(node->ident.tcond);
429     DeleteByteCode(node->ident.fcond);
430     break;
431    
432     case BC_LOCAL:
433    
434     DeleteByteCode(node->local.tcond);
435     DeleteByteCode(node->local.fcond);
436     break;
437     }
438    
439     free(node);
440     }
441     }
442    
443     static void FatalMemoryError(void)
444     {
445     obErrorLog.ThrowError(__FUNCTION__, "Unable to allocate memory for biomolecule residue / chain perception.", obError);
446     // exit(1);
447     }
448    
449     void GenerateByteCodes(ByteCode **node, int resid, int curr, int prev, int bond)
450     {
451     StackType neighbour[4];
452     StackType original;
453     int count,i,j;
454     ByteCode *ptr;
455     bool done,found;
456    
457     if( curr != prev )
458     {
459     if( MonoAtom[curr].atomid < ATOMMINAMINO )
460     {
461     found = false;
462     while( *node && ((*node)->type==BC_IDENT) )
463     {
464     if( (*node)->ident.value == MonoAtom[curr].atomid )
465     {
466     node = (ByteCode**)&(*node)->ident.tcond;
467     found = true;
468     break;
469     }
470     else
471     node = (ByteCode**)&(*node)->ident.fcond;
472     }
473    
474     if (!found)
475     {
476     ptr = AllocateByteCode(BC_IDENT);
477     ptr->ident.tcond = (ByteCode*)0;
478     ptr->ident.fcond = *node;
479     *node = ptr;
480     node = (ByteCode**)&ptr->ident.tcond;
481     ptr->ident.value = MonoAtom[curr].atomid;
482     }
483     MonoBond[bond].index = BondIndex++;
484     done = true;
485     }
486     else if( MonoAtom[curr].index != -1 )
487     {
488     while( *node && ((*node)->type==BC_IDENT) )
489     node = (ByteCode**)&(*node)->ident.fcond;
490    
491     found = false;
492     while( *node && ((*node)->type==BC_LOCAL) )
493     {
494     if( (*node)->local.value == MonoAtom[curr].index )
495     {
496     node = (ByteCode**)&(*node)->local.tcond;
497     found = true;
498     break;
499     }
500     else
501     node = (ByteCode**)&(*node)->local.fcond;
502     }
503    
504     if (!found)
505     {
506     ptr = AllocateByteCode(BC_LOCAL);
507     ptr->local.tcond = (ByteCode*)0;
508     ptr->local.fcond = *node;
509     *node = ptr;
510     node = (ByteCode**)&ptr->local.tcond;
511     ptr->local.value = MonoAtom[curr].index;
512     }
513    
514     MonoBond[bond].index = BondIndex++;
515     done = true;
516     }
517     else
518     {
519     while( *node && ((*node)->type==BC_IDENT) )
520     node = (ByteCode**)&(*node)->ident.fcond;
521     while( *node && ((*node)->type==BC_LOCAL) )
522     node = (ByteCode**)&(*node)->local.fcond;
523    
524     found = false;
525     while( *node && ((*node)->type==BC_ELEM) )
526     {
527     if( (*node)->elem.value == MonoAtom[curr].elem )
528     {
529     node = (ByteCode**)&(*node)->elem.tcond;
530     found = true;
531     break;
532     }
533     else
534     node = (ByteCode**)&(*node)->elem.fcond;
535     }
536    
537     if( !found )
538     {
539     ptr = AllocateByteCode(BC_ELEM);
540     ptr->elem.tcond = (ByteCode*)0;
541     ptr->elem.fcond = *node;
542     *node = ptr;
543     node = (ByteCode**)&ptr->elem.tcond;
544     ptr->elem.value = MonoAtom[curr].elem;
545     }
546    
547     MonoAtom[curr].index = AtomIndex++;
548     MonoBond[bond].index = BondIndex++;
549     done = false;
550     }
551     }
552     else
553     {
554     MonoAtom[curr].index = AtomIndex++;
555     done = false;
556     }
557    
558     count = 0;
559     if (!done)
560     {
561     for( i=0; i<MonoBondCount; i++ )
562     {
563     if( MonoBond[i].src == curr )
564     {
565     if( MonoBond[i].dst != prev )
566     {
567     neighbour[count].atom = MonoBond[i].dst;
568     neighbour[count].bond = i;
569     count++;
570     }
571     }
572     else if( MonoBond[i].dst == curr )
573     {
574     if( MonoBond[i].src != prev )
575     {
576     neighbour[count].atom = MonoBond[i].src;
577     neighbour[count].bond = i;
578     count++;
579     }
580     }
581     }
582    
583     if ( *node && ((*node)->type==BC_EVAL) )
584     {
585     found = false;
586     node = (ByteCode**)&(*node)->eval.next;
587     while( *node && ((*node)->type==BC_COUNT) )
588     {
589     if( (*node)->count.value == count )
590     {
591     node = (ByteCode**)&(*node)->count.tcond;
592     found = true;
593     break;
594     }
595     else
596     node = (ByteCode**)&(*node)->count.fcond;
597     }
598    
599     if( !found )
600     {
601     ptr = AllocateByteCode(BC_COUNT);
602     ptr->count.tcond = (ByteCode*)0;
603     ptr->count.fcond = *node;
604     *node = ptr;
605     node = (ByteCode**)&ptr->count.tcond;
606     ptr->count.value = count;
607     }
608     }
609     else if( count || StrictFlag || StackPtr )
610     {
611     ptr = AllocateByteCode(BC_EVAL);
612     ptr->eval.next = *node;
613     *node = ptr;
614     node = (ByteCode**)&ptr->eval.next;
615    
616     ptr = AllocateByteCode(BC_COUNT);
617     ptr->count.tcond = (ByteCode*)0;
618     ptr->count.fcond = *node;
619     *node = ptr;
620     node = (ByteCode**)&ptr->count.tcond;
621     ptr->count.value = count;
622     }
623     }
624    
625     if( count == 1 )
626     {
627     GenerateByteCodes(node,resid,neighbour[0].atom, curr,neighbour[0].bond);
628     }
629     else if( count == 2 )
630     {
631     original = Stack[StackPtr++];
632     Stack[StackPtr-1] = neighbour[0];
633     Stack[StackPtr-1].prev = curr;
634     GenerateByteCodes(node,resid,neighbour[1].atom,
635     curr,neighbour[1].bond);
636     Stack[StackPtr-1] = neighbour[1];
637     Stack[StackPtr-1].prev = curr;
638     GenerateByteCodes(node,resid,neighbour[0].atom,
639     curr,neighbour[0].bond);
640     Stack[--StackPtr] = original;
641     }
642     else if( count )
643     {
644     #ifdef HAVE_SSTREAM
645     stringstream errorMsg;
646     #else
647     strstream errorMsg;
648     #endif
649     errorMsg << "Maximum Monomer Fanout Exceeded!" << endl;
650     errorMsg << "Residue " << ChainsResName[resid] << " atom "
651     << curr << endl;
652     errorMsg << "Previous = " << prev << " Fanout = " << count << endl;
653     obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obWarning);
654     // exit(1);
655     }
656     else if( StackPtr )
657     {
658     StackPtr--;
659     GenerateByteCodes(node,resid,Stack[StackPtr].atom,
660     Stack[StackPtr].prev,Stack[StackPtr].bond);
661     StackPtr++;
662     }
663     else if( !(*node) )
664     {
665     ptr = AllocateByteCode(BC_ASSIGN);
666     ptr->assign.resid = resid;
667     ptr->assign.atomid = (int*)malloc(AtomIndex*sizeof(int));
668     if( !ptr->assign.atomid )
669     FatalMemoryError();
670     for( i=0; i<MonoAtomCount; i++ )
671     if( (j=MonoAtom[i].index) != -1 )
672     ptr->assign.atomid[j] = MonoAtom[i].atomid;
673     if( BondIndex )
674     {
675     ptr->assign.bflags = (int*)malloc(BondIndex*sizeof(int));
676     for( i=0; i<MonoBondCount; i++ )
677     if( (j=MonoBond[i].index) != -1 )
678     ptr->assign.bflags[j] = MonoBond[i].flag;
679     }
680     *node = ptr;
681     }
682     else if( (*node)->type == BC_ASSIGN )
683     {
684     if( (*node)->assign.resid != resid )
685     {
686     fputs("Error: Duplicated Monomer Specification!\n",stderr);
687     fprintf(stderr,"Residue %s matches resid",ChainsResName[resid]);
688     fprintf(stderr,"ue %s!\n",ChainsResName[(*node)->assign.resid]);
689     }
690     }
691    
692     /* Restore State! */
693     if( curr != prev )
694     {
695     if( !done )
696     {
697     MonoAtom[curr].index = -1;
698     AtomIndex--;
699     }
700     MonoBond[bond].index = -1;
701     BondIndex--;
702     }
703     }
704    
705     //////////////////////////////////////////////////////////////////////////////
706     // Constructors / Destructors
707     //////////////////////////////////////////////////////////////////////////////
708    
709     // validated
710     OBChainsParser::OBChainsParser(void)
711     {
712     int i, res = RESIDMIN;
713    
714     PDecisionTree = (ByteCode*)0;
715     for( i=0 ; i < AMINOMAX ; i++ )
716     {
717     strcpy(ChainsResName[res],AminoAcids[i].name);
718     DefineMonomer(&PDecisionTree,res,AminoAcids[i].data);
719     res++;
720     }
721    
722     NDecisionTree = (ByteCode*)0;
723     for( i=0 ; i< NUCLEOMAX ; i++ )
724     {
725     strcpy(ChainsResName[res],Nucleotides[i].name);
726     DefineMonomer(&NDecisionTree,res,Nucleotides[i].data);
727     res++;
728     }
729    
730     bitmasks = NULL;
731     hetflags = NULL;
732     atomids = NULL;
733     resids = NULL;
734     resnos = NULL;
735     sernos = NULL;
736     hcounts = NULL;
737     chains = NULL;
738     flags = NULL;
739     }
740    
741     OBChainsParser::~OBChainsParser(void)
742     {
743     DeleteByteCode((ByteCode*)PDecisionTree);
744     DeleteByteCode((ByteCode*)NDecisionTree);
745     }
746    
747     //////////////////////////////////////////////////////////////////////////////
748     // Setup / Cleanup Functions
749     //////////////////////////////////////////////////////////////////////////////
750    
751     void OBChainsParser::SetupMol(OBMol &mol)
752     {
753     CleanupMol();
754    
755     int i;
756     int asize = mol.NumAtoms();
757     int bsize = mol.NumBonds();
758    
759     bitmasks = new unsigned short[asize];
760     resids = new unsigned char[asize];
761     flags = new unsigned char[bsize];
762     hetflags = new bool[asize];
763     atomids = new short[asize];
764     resnos = new short[asize];
765     sernos = new short[asize];
766     hcounts = new char[asize];
767     chains = new char[asize];
768    
769     memset(bitmasks, 0, sizeof(unsigned short) * asize);
770     memset(resids, 0, sizeof(unsigned char) * asize);
771     memset(hetflags, 0, sizeof(bool) * asize);
772     memset(resnos, 0, sizeof(short) * asize);
773     memset(sernos, 0, sizeof(short) * asize);
774     memset(hcounts, 0, sizeof(char) * asize);
775     memset(chains, ' ', sizeof(char) * asize);
776    
777     memset(flags, 0, sizeof(unsigned char) * bsize);
778    
779     for ( i = 0 ; i < asize ; i++ )
780     {
781     atomids[i] = -1;
782     }
783     }
784    
785     void OBChainsParser::CleanupMol(void)
786     {
787     if (bitmasks != NULL)
788     {
789     delete bitmasks;
790     bitmasks = NULL;
791     }
792     if (hetflags != NULL)
793     {
794     delete hetflags;
795     hetflags = NULL;
796     }
797     if (atomids != NULL)
798     {
799     delete atomids;
800     atomids = NULL;
801     }
802     if (resids != NULL)
803     {
804     delete resids;
805     resids = NULL;
806     }
807     if (resnos != NULL)
808     {
809     delete resnos;
810     resnos = NULL;
811     }
812     if (sernos != NULL)
813     {
814     delete sernos;
815     sernos = NULL;
816     }
817     if (hcounts != NULL)
818     {
819     delete hcounts;
820     hcounts = NULL;
821     }
822     if (chains != NULL)
823     {
824     delete chains;
825     chains = NULL;
826     }
827     if (flags != NULL)
828     {
829     delete flags;
830     flags = NULL;
831     }
832     }
833    
834     void OBChainsParser::ClearResidueInformation(OBMol &mol)
835     {
836     OBResidue *residue;
837     vector<OBResidue*> residues;
838     vector<OBResidue*>::iterator r;
839    
840     for (residue = mol.BeginResidue(r) ; residue ; residue = mol.NextResidue(r))
841     residues.push_back(residue);
842    
843     for ( unsigned int i = 0 ; i < residues.size() ; i++ )
844     mol.DeleteResidue(residues[i]);
845    
846     residues.clear();
847     }
848    
849     void OBChainsParser::SetResidueInformation(OBMol &mol, bool nukeSingleResidue)
850     {
851     char buffer[256];
852     char *symbol;
853     string atomid, name;
854    
855     OBAtom *atom;
856     OBResidue *residue;
857     map<short, OBResidue *> resmap;
858    
859     int size = mol.NumAtoms();
860     for ( int i = 0 ; i < size ; i++ )
861     {
862     atom = mol.GetAtom(i+1); // WARNING: ATOM INDEX ISSUE
863    
864     if (atomids[i] == -1)
865     {
866     symbol = etab.GetSymbol(atom->GetAtomicNum());
867     if( symbol[1] )
868     {
869     buffer[0] = symbol[0];
870     buffer[1] = (char) toupper(symbol[1]);
871     }
872     else
873     {
874     buffer[0] = ' ';
875     buffer[1] = symbol[0];
876     }
877     buffer[2] = ' ';
878     buffer[3] = ' ';
879     buffer[4] = '\0';
880     }
881     else if (atom->IsHydrogen())
882     {
883     if (hcounts[i])
884     sprintf(buffer, "%cH%.2s", hcounts[i]+'0', ChainsAtomName[atomids[i]]+2);
885     else
886     sprintf(buffer, "H%.2s", ChainsAtomName[atomids[i]]+2);
887     }
888     else
889     sprintf(buffer, "%.4s", ChainsAtomName[atomids[i]]);
890    
891     if (buffer[3] == ' ')
892     buffer[3] = '\0';
893    
894     atomid = (buffer[0] == ' ') ? buffer + 1 : buffer;
895    
896     if (resmap.find(resnos[i]) != resmap.end())
897     {
898     residue = resmap[resnos[i]];
899     residue->AddAtom(atom);
900     residue->SetAtomID(atom, atomid);
901     residue->SetHetAtom(atom, hetflags[i]);
902     residue->SetSerialNum(atom, sernos[i]);
903     }
904     else
905     {
906     name = ChainsResName[resids[i]];
907     residue = mol.NewResidue();
908    
909     residue->SetName(name);
910     residue->SetNum(resnos[i]);
911     residue->SetChain(chains[i]);
912     residue->SetChainNum((chains[i] > 'A') ? (int)(chains[i] - 'A') : 1);
913    
914     residue->AddAtom(atom);
915     residue->SetAtomID(atom, atomid);
916     residue->SetHetAtom(atom, hetflags[i]);
917     residue->SetSerialNum(atom, sernos[i]);
918    
919     resmap[resnos[i]] = residue;
920     }
921     }
922    
923     if (mol.NumResidues() == 1 && nukeSingleResidue)
924     mol.DeleteResidue(mol.GetResidue(0));
925     }
926    
927     ////////////////////////////////////////////////////////////////////////////////
928     // Perception Functions
929     ////////////////////////////////////////////////////////////////////////////////
930    
931     bool OBChainsParser::PerceiveChains(OBMol &mol, bool nukeSingleResidue)
932     {
933     bool result = true;
934    
935     SetupMol(mol);
936     ClearResidueInformation(mol);
937    
938     result = DetermineHetAtoms(mol) && result;
939     result = DetermineConnectedChains(mol) && result;
940     result = DeterminePeptideBackbone(mol) && result;
941     result = DeterminePeptideSidechains(mol) && result;
942     result = DetermineNucleicBackbone(mol) && result;
943     result = DetermineNucleicSidechains(mol) && result;
944     result = DetermineHydrogens(mol) && result;
945    
946     SetResidueInformation(mol, nukeSingleResidue);
947     CleanupMol();
948    
949     obErrorLog.ThrowError(__FUNCTION__,
950     "Ran OpenBabel::PerceiveChains", obAuditMsg);
951    
952     return result;
953     }
954    
955     ////////////////////////////////////////////////////////////////////////////////
956     // Hetero Atom Perception
957     ////////////////////////////////////////////////////////////////////////////////
958    
959     bool OBChainsParser::DetermineHetAtoms(OBMol &mol)
960     {
961     OBAtom *atom;
962     vector<OBNodeBase *>::iterator a;
963     for (atom = mol.BeginAtom(a) ; atom ; atom = mol.NextAtom(a))
964     if (!atom->IsHydrogen() && atom->GetValence() == 0)
965     {
966     resids[atom->GetIdx()-1] = (atom->IsOxygen()) ? 1 : 2;
967     hetflags[atom->GetIdx()-1] = true;
968     }
969     return true;
970     }
971    
972     ////////////////////////////////////////////////////////////////////////////////
973     // Connected Chain Perception
974     ////////////////////////////////////////////////////////////////////////////////
975    
976     bool OBChainsParser::DetermineConnectedChains(OBMol &mol)
977     {
978     int resid;
979     int resno;
980     int count;
981     int size;
982     int i,idx;
983     int numAtoms;
984    
985     resno = 1;
986     count = 0;
987     numAtoms = mol.NumAtoms();
988    
989     OBAtom *atom;
990     vector<OBNodeBase *>::iterator a;
991     for (atom = mol.BeginAtom(a) ; atom ; atom = mol.NextAtom(a))
992     {
993     idx = atom->GetIdx() - 1;
994     if (!hetflags[idx] && chains[idx] == ' ' && !atom->IsHydrogen())
995     {
996     size = RecurseChain(mol, idx, 'A' + count);
997     if (size < 10)
998     {
999     if (size == 1 && atom->IsOxygen())
1000     resid = 1; /* HOH */
1001     else
1002     resid = 2;
1003    
1004     for (i = 0 ; i < numAtoms ; i++)
1005     {
1006     if (chains[i] == ('A' + count))
1007     {
1008     hetflags[i] = true;
1009     resids[i] = resid;
1010     resnos[i] = resno;
1011     chains[i] = ' ';
1012     }
1013     }
1014     resno++;
1015     }
1016     else
1017     count++;
1018     }
1019     }
1020    
1021     if( count == 1 )
1022     for ( i = 0 ; i < numAtoms ; i++ )
1023     chains[i] = ' ';
1024    
1025     return true;
1026     }
1027    
1028     int OBChainsParser::RecurseChain(OBMol &mol, int i, int c)
1029     {
1030     OBAtom *atom, *nbr;
1031     vector<OBEdgeBase *>::iterator b;
1032     int result, index;
1033    
1034     atom = mol.GetAtom(i+1);
1035     if (atom->IsHydrogen() )
1036     return 0;
1037    
1038     result = 1;
1039     chains[i] = c;
1040    
1041     for (nbr = atom->BeginNbrAtom(b) ; nbr ; nbr = atom->NextNbrAtom(b))
1042     {
1043     index = nbr->GetIdx() - 1;
1044     if (chains[index] == ' ')
1045     result += RecurseChain(mol, index,c);
1046     }
1047    
1048     return (result);
1049     }
1050    
1051     ////////////////////////////////////////////////////////////////////////////////
1052     // Peptide Backbone Perception
1053     ////////////////////////////////////////////////////////////////////////////////
1054    
1055     bool OBChainsParser::DeterminePeptideBackbone(OBMol &mol)
1056     {
1057     ConstrainBackbone(mol, Peptide, MAXPEPTIDE);
1058    
1059     int i, max = mol.NumAtoms();
1060    
1061     /*
1062     int count = 0;
1063     for ( i = 0 ; i < max ; i++ )
1064     if ( bitmasks[i]&BitCAAll )
1065     count++;
1066    
1067     fprintf(stderr,"%d alpha carbons\n",count);
1068     */
1069    
1070     /* Order Peptide Backbone */
1071    
1072     for ( i = 0 ; i < max ; i++ )
1073     if (atomids[i] == -1)
1074     {
1075     if( bitmasks[i] & BitNTer )
1076     {
1077     atomids[i] = AI_N;
1078     TracePeptideChain(mol,i,1);
1079     }
1080     else if( (bitmasks[i]&BitNPT) && !(bitmasks[i]&BitN) )
1081     {
1082     atomids[i] = AI_N;
1083     TracePeptideChain(mol,i,1);
1084     }
1085     }
1086    
1087     /* Carbonyl Double Bond */
1088    
1089     OBBond *bond;
1090     vector<OBEdgeBase*>::iterator b;
1091     for (bond = mol.BeginBond(b) ; bond ; bond = mol.NextBond(b))
1092     {
1093     if ((atomids[bond->GetBeginAtomIdx()-1] == 2 && atomids[bond->GetEndAtomIdx()-1] == 3) ||
1094     (atomids[bond->GetBeginAtomIdx()-1] == 3 && atomids[bond->GetEndAtomIdx()-1] == 2))
1095     flags[bond->GetIdx()] |= BF_DOUBLE;
1096     }
1097    
1098     return true;
1099     }
1100    
1101     void OBChainsParser::ConstrainBackbone(OBMol &mol, Template *templ, int tmax)
1102     {
1103     static OBAtom *neighbour[6];
1104     Template *pep;
1105     OBAtom *na,*nb,*nc,*nd;
1106     OBAtom *atom, *nbr;
1107     bool change, result;
1108     int count;
1109     int i,idx;
1110    
1111     vector<OBNodeBase *>::iterator a;
1112     vector<OBEdgeBase *>::iterator b;
1113    
1114     /* First Pass */
1115    
1116     for (atom = mol.BeginAtom(a) ; atom ; atom = mol.NextAtom(a))
1117     {
1118     idx = atom->GetIdx() - 1;
1119     bitmasks[idx] = 0;
1120     for ( i = 0 ; i < tmax ; i++ )
1121     if ( (static_cast<unsigned int>(templ[i].elem) == atom->GetAtomicNum())
1122     &&
1123     (static_cast<unsigned int>(templ[i].count) == atom->GetValence()))
1124     bitmasks[idx] |= templ[i].flag;
1125     }
1126    
1127     /* Second Pass */
1128    
1129     do
1130     {
1131     change = false;
1132     for (atom = mol.BeginAtom(a) ; atom ; atom = mol.NextAtom(a))
1133     {
1134     idx = atom->GetIdx() - 1;
1135     if (bitmasks[idx]) // Determine Neighbours
1136     {
1137     count = 0;
1138     for (nbr = atom->BeginNbrAtom(b) ; nbr ; nbr = atom->NextNbrAtom(b))
1139     if (!nbr->IsHydrogen())
1140     neighbour[count++] = nbr;
1141    
1142     na = neighbour[0];
1143     nb = neighbour[1];
1144     nc = neighbour[2];
1145     nd = neighbour[3];
1146    
1147     for ( i = 0 ; i < tmax ; i++ )
1148     if ( templ[i].flag & bitmasks[idx] )
1149     {
1150     pep = &templ[i];
1151     result = true;
1152    
1153     if (count == 4)
1154     result = Match4Constraints(pep,na,nb,nc,nd);
1155     else if (count == 3)
1156     result = Match3Constraints(pep,na,nb,nc);
1157     else if (count == 2)
1158     result = Match2Constraints(pep,na,nb);
1159     else // count == 1
1160     result = MatchConstraint(na,pep->n1);
1161    
1162     if(result == false)
1163     {
1164     bitmasks[idx] &= ~pep->flag;
1165     change = true;
1166     }
1167     }
1168     }
1169     }
1170     }
1171     while( change );
1172     }
1173    
1174     bool OBChainsParser::MatchConstraint(OBAtom *atom, int mask)
1175     {
1176     if (atom == NULL)
1177     return (false);
1178    
1179     if( mask < 0 )
1180     return(atom->GetAtomicNum() == static_cast<unsigned int>(-mask));
1181     else
1182     return(((bitmasks[atom->GetIdx()-1]&mask) == 0) ? false : true);
1183     }
1184    
1185     bool OBChainsParser::Match2Constraints(Template *tmpl, OBAtom *na, OBAtom *nb)
1186     {
1187     if (na == NULL || nb == NULL)
1188     return (false); // don't even try to evaluate it
1189    
1190     if( MatchConstraint(na,tmpl->n2) )
1191     if( MatchConstraint(nb,tmpl->n1) )
1192     return( true );
1193     if( MatchConstraint(nb,tmpl->n2) )
1194     if( MatchConstraint(na,tmpl->n1) )
1195     return( true );
1196     return( false );
1197     }
1198    
1199     bool OBChainsParser::Match3Constraints(Template *tmpl, OBAtom *na, OBAtom *nb, OBAtom *nc)
1200     {
1201     if (na == NULL || nb == NULL || nc == NULL)
1202     return (false); // don't even try to evaluate it
1203    
1204     if( MatchConstraint(na,tmpl->n3) )
1205     if( Match2Constraints(tmpl,nb,nc) )
1206     return( true );
1207     if( MatchConstraint(nb,tmpl->n3) )
1208     if( Match2Constraints(tmpl,na,nc) )
1209     return( true );
1210     if( MatchConstraint(nc,tmpl->n3) )
1211     if( Match2Constraints(tmpl,na,nb) )
1212     return( true );
1213     return( false );
1214     }
1215    
1216     bool OBChainsParser::Match4Constraints(Template *tmpl, OBAtom *na, OBAtom *nb, OBAtom *nc, OBAtom *nd)
1217     {
1218     if (na == NULL || nb == NULL || nc == NULL || nd == NULL)
1219     return (false); // don't even try to evaluate it
1220    
1221     if( MatchConstraint(na,tmpl->n4) )
1222     if( Match3Constraints(tmpl,nb,nc,nd) )
1223     return( true );
1224     if( MatchConstraint(nb,tmpl->n4) )
1225     if( Match3Constraints(tmpl,na,nc,nd) )
1226     return( true );
1227     if( MatchConstraint(nc,tmpl->n4) )
1228     if( Match3Constraints(tmpl,na,nb,nd) )
1229     return( true );
1230     if( MatchConstraint(nd,tmpl->n4) )
1231     if( Match3Constraints(tmpl,na,nb,nc) )
1232     return( true );
1233     return( false );
1234     }
1235    
1236     void OBChainsParser::TracePeptideChain(OBMol &mol, int i, int r)
1237     {
1238     int neighbour[4];
1239     int na,nb,nc;
1240     OBAtom *atom, *nbr;
1241     int count;
1242     int j,k,idx;
1243    
1244     vector<OBEdgeBase *>::iterator b;
1245    
1246     /* Determine Neighbours */
1247    
1248     atom = mol.GetAtom(i+1);
1249     idx = atom->GetIdx() - 1;
1250    
1251     count = 0;
1252     for (nbr = atom->BeginNbrAtom(b) ; nbr ; nbr = atom->NextNbrAtom(b))
1253     if (!nbr->IsHydrogen())
1254     neighbour[count++] = nbr->GetIdx()-1;
1255    
1256     resnos[idx] = r;
1257    
1258     na = neighbour[0];
1259     nb = neighbour[1];
1260     nc = neighbour[2];
1261    
1262     switch( atomids[i] )
1263     {
1264     case(AI_N):
1265     for( j=0; j<count; j++ )
1266     if( bitmasks[neighbour[j]] & BitCAAll )
1267     {
1268     atomids[neighbour[j]] = AI_CA;
1269     TracePeptideChain(mol,neighbour[j],r);
1270     }
1271     break;
1272    
1273     case(AI_CA):
1274     if( count == 3 )
1275     {
1276     if ( bitmasks[na] & BitNAll )
1277     na = nc;
1278     else if ( bitmasks[nb] & BitNAll )
1279     nb = nc;
1280    
1281     if ( bitmasks[na] & BitC )
1282     {
1283     j = na;
1284     k = nb;
1285     }
1286     else if ( bitmasks[nb] & BitC )
1287     {
1288     j = nb;
1289     k = na;
1290     }
1291     else if( bitmasks[na] & BitCAll )
1292     {
1293     j = na;
1294     k = nb;
1295     }
1296     else /* bitmasks[nb] & BitCAll */
1297     {
1298     j = nb;
1299     k = na;
1300     }
1301    
1302     atomids[j] = AI_C;
1303     bitmasks[k] = 0;
1304    
1305     TracePeptideChain(mol,j,r);
1306     }
1307     else /* count == 2 */
1308     {
1309     if ( bitmasks[na] & BitCAll )
1310     {
1311     atomids[na] = AI_C;
1312     TracePeptideChain(mol,na,r);
1313     }
1314     else
1315     {
1316     atomids[nb] = AI_C;
1317     TracePeptideChain(mol,nb,r);
1318     }
1319     }
1320     break;
1321    
1322     case(AI_C):
1323     k = AI_O;
1324     for ( j = 0; j < count; j++ )
1325     {
1326     if ( bitmasks[neighbour[j]] & BitNAll )
1327     {
1328     atomids[neighbour[j]] = AI_N;
1329     TracePeptideChain(mol,neighbour[j],r+1);
1330     }
1331     else if( bitmasks[neighbour[j]] & BitOAll )
1332     {
1333     atomids[neighbour[j]] = k;
1334     resnos[neighbour[j]] = r;
1335     k = AI_OXT; /* OXT */
1336     }
1337     }
1338     break;
1339     }
1340     }
1341    
1342     ////////////////////////////////////////////////////////////////////////////////
1343     // Peptide Sidechains Perception
1344     ////////////////////////////////////////////////////////////////////////////////
1345    
1346     bool OBChainsParser::DeterminePeptideSidechains(OBMol &mol)
1347     {
1348     int resid;
1349     int max = mol.NumAtoms();
1350    
1351     for (int i = 0 ; i < max ; i++)
1352     if (atomids[i] == 1)
1353     {
1354     resid = IdentifyResidue(PDecisionTree, mol, i, resnos[i]);
1355     AssignResidue(mol,resnos[i],chains[i],resid);
1356     }
1357    
1358     return true;
1359     }
1360    
1361     void OBChainsParser::AssignResidue(OBMol &mol, int r, int c, int i)
1362     {
1363     int max = mol.NumAtoms();
1364     for (int j = 0 ; j < max ; j++)
1365     if ((resnos[j] == r) && (chains[j] == c) && !hetflags[j])
1366     resids[j] = i;
1367     }
1368    
1369     int OBChainsParser::IdentifyResidue(void *tree, OBMol &mol, int seed, int resno)
1370     {
1371     ByteCode *ptr;
1372    
1373     int AtomCount, BondCount;
1374     int curr,prev,bond;
1375     int bcount;
1376     int i,j;
1377    
1378     ptr = (ByteCode *) tree;
1379     bcount = 0;
1380    
1381     Stack[0].atom = seed;
1382     Stack[0].prev = seed;
1383     StackPtr = 0;
1384    
1385     ResMonoAtom[0] = seed;
1386     AtomCount = 1;
1387     BondCount = 0;
1388    
1389     OBAtom *atom, *nbr;
1390     vector<OBEdgeBase *>::iterator b;
1391    
1392     while( ptr )
1393     switch(ptr->type)
1394     {
1395     case(BC_IDENT): curr = Stack[StackPtr-1].atom;
1396     if( atomids[curr] == ptr->ident.value )
1397     {
1398     bond = Stack[StackPtr-1].bond;
1399     ResMonoBond[BondCount++] = bond;
1400     ptr = ptr->ident.tcond;
1401     StackPtr--;
1402     }
1403     else
1404     ptr = ptr->ident.fcond;
1405     break;
1406    
1407     case(BC_LOCAL): curr = Stack[StackPtr-1].atom;
1408     if( curr == ResMonoAtom[ptr->local.value] )
1409     {
1410     bond = Stack[StackPtr-1].bond;
1411     ResMonoBond[BondCount++] = bond;
1412     ptr = ptr->local.tcond;
1413     StackPtr--;
1414     }
1415     else
1416     ptr = ptr->local.fcond;
1417     break;
1418    
1419     case(BC_ELEM): curr = Stack[StackPtr-1].atom;
1420     if( mol.GetAtom(curr+1)->GetAtomicNum() == static_cast<unsigned int>(ptr->elem.value)
1421     )
1422     {
1423     bond = Stack[StackPtr-1].bond;
1424     ResMonoAtom[AtomCount++] = curr;
1425     ResMonoBond[BondCount++] = bond;
1426     resnos[curr] = resno;
1427     ptr = ptr->elem.tcond;
1428     StackPtr--;
1429     }
1430     else
1431     ptr = ptr->elem.fcond;
1432     break;
1433    
1434     case(BC_EVAL): bcount = 0;
1435     curr = Stack[StackPtr].atom;
1436     prev = Stack[StackPtr].prev;
1437    
1438     atom = mol.GetAtom(curr+1);
1439     for (nbr = atom->BeginNbrAtom(b); nbr; nbr = atom->NextNbrAtom(b))
1440     {
1441     j = nbr->GetIdx() - 1;
1442     if (!((curr == prev) && bitmasks[j]) && (j != prev) && !(nbr->IsHydrogen()))
1443     {
1444     Stack[StackPtr].prev = curr;
1445     Stack[StackPtr].atom = j;
1446     Stack[StackPtr].bond = (*b)->GetIdx();
1447     StackPtr++;
1448     bcount++;
1449     }
1450     }
1451    
1452     ptr = ptr->eval.next;
1453     break;
1454    
1455     case(BC_COUNT):
1456     if( bcount == ptr->count.value )
1457     {
1458     ptr = ptr->count.tcond;
1459     }
1460     else
1461     ptr = ptr->count.fcond;
1462     break;
1463    
1464     case(BC_ASSIGN):
1465     for( i=0; i<AtomCount; i++ )
1466     if( !bitmasks[ResMonoAtom[i]] )
1467     {
1468     j = ptr->assign.atomid[i];
1469     atomids[ResMonoAtom[i]] = j;
1470     }
1471     for( i=0; i<BondCount; i++ )
1472     {
1473     j = ptr->assign.bflags[i];
1474     flags[ResMonoBond[i]] = j;
1475     }
1476     return( ptr->assign.resid );
1477     break;
1478    
1479     default: /* Illegal Instruction! */
1480     return( 0 );
1481     }
1482     return 0;
1483     }
1484    
1485     ////////////////////////////////////////////////////////////////////////////////
1486     // Nucleic Backbone Perception
1487     ////////////////////////////////////////////////////////////////////////////////
1488    
1489     bool OBChainsParser::DetermineNucleicBackbone(OBMol &mol)
1490     {
1491     ConstrainBackbone(mol, Nucleotide, MAXNUCLEIC);
1492    
1493     int i, max = mol.NumAtoms();
1494    
1495     /*
1496     int count = 0;
1497     for ( i = 0 ; i < max ; i++ )
1498     if ( bitmasks[i] & BitC5 )
1499     count++;
1500    
1501     fprintf(stderr,"%d sugar phosphates\n",count);
1502     */
1503    
1504     /* Order Nucleic Backbone */
1505    
1506     for( i = 0 ; i < max ; i++ )
1507     if( atomids[i] == -1 )
1508     {
1509     if( bitmasks[i] & BitPTer )
1510     {
1511     atomids[i] = AI_P;
1512     TraceNucleicChain(mol,i,1);
1513     }
1514     else if( bitmasks[i] & BitO5Ter )
1515     {
1516     atomids[i] = AI_O5;
1517     TraceNucleicChain(mol,i,1);
1518     }
1519     }
1520    
1521     return true;
1522     }
1523    
1524     void OBChainsParser::TraceNucleicChain(OBMol &mol, int i, int r)
1525     {
1526     int neighbour[4];
1527     int na,nb,nc;
1528     int count;
1529     int j,k;
1530    
1531     OBAtom *atom, *nbr;
1532     vector<OBEdgeBase *>::iterator b;
1533    
1534     count = 0;
1535     atom = mol.GetAtom(i + 1);
1536     for (nbr = atom->BeginNbrAtom(b) ; nbr ; nbr = atom->NextNbrAtom(b))
1537     if (!nbr->IsHydrogen())
1538     neighbour[count++] = nbr->GetIdx() - 1;
1539    
1540     resnos[i] = r;
1541    
1542     na = neighbour[0];
1543     nb = neighbour[1];
1544     nc = neighbour[2];
1545    
1546     switch( atomids[i] )
1547     {
1548     case(AI_P):
1549     k = AI_O1P; /* O1P */
1550     for( j=0; j<count; j++ )
1551     {
1552     if( bitmasks[neighbour[j]] & BitO5 )
1553     {
1554     atomids[neighbour[j]] = AI_O5;
1555     TraceNucleicChain(mol,neighbour[j],r);
1556     }
1557     else if( bitmasks[neighbour[j]] & BitOP )
1558     {
1559     atomids[neighbour[j]] = k;
1560     resnos[neighbour[j]] = r;
1561     k = AI_O2P; /* O2P */
1562     }
1563     }
1564    
1565     break;
1566    
1567     case(AI_O5):
1568     for( j=0; j<count; j++ )
1569     if( bitmasks[neighbour[j]] & BitC5 )
1570     {
1571     atomids[neighbour[j]] = AI_C5;
1572     TraceNucleicChain(mol,neighbour[j],r);
1573     }
1574    
1575     break;
1576    
1577     case(AI_C5):
1578     for( j=0 ; j<count; j++ )
1579     if( bitmasks[neighbour[j]] & BitC4 )
1580     {
1581     atomids[neighbour[j]] = AI_C4;
1582     TraceNucleicChain(mol,neighbour[j],r);
1583     }
1584    
1585     break;
1586    
1587     case(AI_C4):
1588     for( j=0; j<count; j++ )
1589     {
1590     if( bitmasks[neighbour[j]] & BitC3 )
1591     {
1592     atomids[neighbour[j]] = AI_C3;
1593     TraceNucleicChain(mol,neighbour[j],r);
1594     }
1595     else if( bitmasks[neighbour[j]] & BitO4 )
1596     {
1597     atomids[neighbour[j]] = AI_O4;
1598     resnos[neighbour[j]] = r;
1599     }
1600     }
1601    
1602     break;
1603    
1604     case(AI_C3):
1605     for( j=0; j<count; j++ )
1606     {
1607     if( bitmasks[neighbour[j]] & BitO3All )
1608     {
1609     atomids[neighbour[j]] = AI_O3;
1610     TraceNucleicChain(mol,neighbour[j],r);
1611     }
1612     else if( bitmasks[neighbour[j]] & BitC2All )
1613     {
1614     atomids[neighbour[j]] = AI_C2;
1615     TraceNucleicChain(mol,neighbour[j],r);
1616     }
1617     }
1618    
1619     break;
1620    
1621     case(AI_O3):
1622     for( j=0; j<count; j++ )
1623     if( bitmasks[neighbour[j]] & BitP )
1624     {
1625     atomids[neighbour[j]] = AI_P;
1626     TraceNucleicChain(mol,neighbour[j],r+1);
1627     }
1628    
1629     break;
1630    
1631     case(AI_C2):
1632     for( j=0; j<count; j++ )
1633     {
1634     if( bitmasks[neighbour[j]] & BitC1 )
1635     {
1636     atomids[neighbour[j]] = AI_C1;
1637     resnos[neighbour[j]] = r;
1638     }
1639     else if( bitmasks[neighbour[j]] & BitO2 )
1640     {
1641     atomids[neighbour[j]] = AI_O2;
1642     resnos[neighbour[j]] = r;
1643     }
1644     }
1645    
1646     break;
1647     }
1648     }
1649    
1650     ////////////////////////////////////////////////////////////////////////////////
1651     // Nucleic Sidechains Perception
1652     ////////////////////////////////////////////////////////////////////////////////
1653    
1654     bool OBChainsParser::DetermineNucleicSidechains(OBMol &mol)
1655     {
1656     for( unsigned int i = 0 ; i < mol.NumAtoms() ; i++ )
1657     if( atomids[i] == 49 )
1658     {
1659     int resid = IdentifyResidue(NDecisionTree,mol,i,resnos[i]);
1660     AssignResidue(mol,resnos[i],chains[i],resid);
1661     }
1662    
1663     return true;
1664     }
1665    
1666     ////////////////////////////////////////////////////////////////////////////////
1667     // Hydrogens Perception
1668     ////////////////////////////////////////////////////////////////////////////////
1669    
1670     bool OBChainsParser::DetermineHydrogens(OBMol &mol)
1671     {
1672     OBAtom *atom, *nbr;
1673     int idx,sidx;
1674    
1675     int max = mol.NumAtoms();
1676     for ( int i = 0 ; i < max ; i++ )
1677     hcounts[i] = 0;
1678    
1679     /* First Pass */
1680    
1681     vector<OBNodeBase*>::iterator a;
1682     vector<OBEdgeBase*>::iterator b;
1683    
1684     for(atom = mol.BeginAtom(a); atom ; atom = mol.NextAtom(a))
1685     if(atom->IsHydrogen())
1686     {
1687     nbr = atom->BeginNbrAtom(b);
1688     if (nbr != NULL)
1689     {
1690     idx = atom->GetIdx() - 1;
1691     sidx = nbr->GetIdx() - 1;
1692    
1693     hcounts[idx] = ++hcounts[sidx];
1694     hetflags[idx] = hetflags[sidx];
1695     atomids[idx] = atomids[sidx];
1696     resids[idx] = resids[sidx];
1697     resnos[idx] = resnos[sidx];
1698     }
1699     }
1700    
1701     /* Second Pass */
1702    
1703     for(atom = mol.BeginAtom(a) ; atom ; atom = mol.NextAtom(a))
1704     if (atom->IsHydrogen())
1705     {
1706     nbr = atom->BeginNbrAtom(b);
1707     if (nbr != NULL && hcounts[nbr->GetIdx()-1] == 1)
1708     hcounts[atom->GetIdx()-1] = 0;
1709     }
1710    
1711     return true;
1712     }
1713    
1714     ////////////////////////////////////////////////////////////////////////////////
1715     // Utility Functions
1716     ////////////////////////////////////////////////////////////////////////////////
1717    
1718     // validated
1719     void OBChainsParser::DefineMonomer(void **tree, int resid, char *smiles)
1720     {
1721     int i;
1722    
1723     MonoAtomCount = 0;
1724     MonoBondCount = 0;
1725    
1726     ParseSmiles(smiles,-1);
1727    
1728     for( i=0; i<MonoBondCount; i++ )
1729     MonoBond[i].index = -1;
1730     for( i=0; i<MonoAtomCount; i++ )
1731     MonoAtom[i].index = -1;
1732     AtomIndex = BondIndex = 0;
1733    
1734     StackPtr = 0;
1735     GenerateByteCodes((ByteCode**)tree, resid, 0, 0, 0 );
1736     }
1737    
1738     int OBChainsParser::IdentifyElement(char *ptr)
1739     {
1740     int ch;
1741    
1742     ch = toupper(ptr[1]);
1743     switch( toupper(ptr[0]) )
1744     {
1745     case(' '): switch( ch )
1746     {
1747     case('B'): return( 5 );
1748     case('C'): return( 6 );
1749     case('D'): return( 1 );
1750     case('F'): return( 9 );
1751     case('H'): return( 1 );
1752     case('I'): return( 53 );
1753     case('K'): return( 19 );
1754     case('L'): return( 1 );
1755     case('N'): return( 7 );
1756     case('O'): return( 8 );
1757     case('P'): return( 15 );
1758     case('S'): return( 16 );
1759     case('U'): return( 92 );
1760     case('V'): return( 23 );
1761     case('W'): return( 74 );
1762     case('Y'): return( 39 );
1763     }
1764     break;
1765    
1766     case('A'): switch( ch )
1767     {
1768     case('C'): return( 89 );
1769     case('G'): return( 47 );
1770     case('L'): return( 13 );
1771     case('M'): return( 95 );
1772     case('R'): return( 18 );
1773     case('S'): return( 33 );
1774     case('T'): return( 85 );
1775     case('U'): return( 79 );
1776     }
1777     break;
1778    
1779     case('B'): switch( ch )
1780     {
1781     case('A'): return( 56 );
1782     case('E'): return( 4 );
1783     case('I'): return( 83 );
1784     case('K'): return( 97 );
1785     case('R'): return( 35 );
1786     case(' '): return( 5 );
1787     }
1788     break;
1789    
1790     case('C'): switch( ch )
1791     {
1792     case('A'): return( 20 );
1793     case('D'): return( 48 );
1794     case('E'): return( 58 );
1795     case('F'): return( 98 );
1796     case('L'): return( 17 );
1797     case('M'): return( 96 );
1798     case('O'): return( 27 );
1799     case('R'): return( 24 );
1800     case('S'): return( 55 );
1801     case('U'): return( 29 );
1802     case(' '): return( 6 );
1803     }
1804     break;
1805    
1806     case('D'): if( ch=='Y' )
1807     {
1808     return( 66 );
1809     }
1810     else if( ch==' ' )
1811     return( 1 );
1812     break;
1813    
1814     case('E'): if( ch=='R' )
1815     {
1816     return( 68 );
1817     }
1818     else if( ch=='S' )
1819     {
1820     return( 99 );
1821     }
1822     else if( ch=='U' )
1823     return( 63 );
1824     break;
1825    
1826     case('F'): if( ch=='E' )
1827     {
1828     return( 26 );
1829     }
1830     else if( ch=='M' )
1831     {
1832     return( 100 );
1833     }
1834     else if( ch=='R' )
1835     {
1836     return( 87 );
1837     }
1838     else if( ch=='F' )
1839     return( 9 );
1840     break;
1841    
1842     case('G'): if( ch=='A' )
1843     {
1844     return( 31 );
1845     }
1846     else if( ch=='D' )
1847     {
1848     return( 64 );
1849     }
1850     else if( ch=='E' )
1851     return( 32 );
1852     break;
1853    
1854     case('H'): if( ch=='E' )
1855     {
1856     return( 2 );
1857     }
1858     else if( ch=='F' )
1859     {
1860     return( 72 );
1861     }
1862     else if( ch=='G' )
1863     {
1864     return( 80 );
1865     }
1866     else if( ch=='O' )
1867     {
1868     return( 67 );
1869     }
1870     else if( ch==' ' )
1871     return( 1 );
1872     break;
1873    
1874     case('I'): if( ch=='N' )
1875     {
1876     return( 49 );
1877     }
1878     else if( ch=='R' )
1879     {
1880     return( 77 );
1881     }
1882     else if( ch==' ' )
1883     return( 53 );
1884     break;
1885    
1886     case('K'): if( ch=='R' )
1887     {
1888     return( 36 );
1889     }
1890     else if( ch==' ' )
1891     return( 19 );
1892     break;
1893    
1894     case('L'): if( ch=='A' )
1895     {
1896     return( 57 );
1897     }
1898     else if( ch=='I' )
1899     {
1900     return( 3 );
1901     }
1902     else if( (ch=='R') || (ch=='W') )
1903     {
1904     return( 103 );
1905     }
1906     else if( ch=='U' )
1907     {
1908     return( 71 );
1909     }
1910     else if( ch==' ' )
1911     return( 1 );
1912     break;
1913    
1914     case('M'): if( ch=='D' )
1915     {
1916     return( 101 );
1917     }
1918     else if( ch=='G' )
1919     {
1920     return( 12 );
1921     }
1922     else if( ch=='N' )
1923     {
1924     return( 25 );
1925     }
1926     else if( ch=='O' )
1927     return( 42 );
1928     break;
1929    
1930     case('N'): switch( ch )
1931     {
1932     case('A'): return( 11 );
1933     case('B'): return( 41 );
1934     case('D'): return( 60 );
1935     case('E'): return( 10 );
1936     case('I'): return( 28 );
1937     case('O'): return( 102 );
1938     case('P'): return( 93 );
1939     case(' '): return( 7 );
1940     }
1941     break;
1942    
1943     case('O'): if( ch=='S' )
1944     {
1945     return( 76 );
1946     }
1947     else if( ch==' ' )
1948     return( 8 );
1949     break;
1950    
1951     case('P'): switch( ch )
1952     {
1953     case('A'): return( 91 );
1954     case('B'): return( 82 );
1955     case('D'): return( 46 );
1956     case('M'): return( 61 );
1957     case('O'): return( 84 );
1958     case('R'): return( 59 );
1959     case('T'): return( 78 );
1960     case('U'): return( 94 );
1961     case(' '): return( 15 );
1962     }
1963     break;
1964    
1965     case('R'): switch( ch )
1966     {
1967     case('A'): return( 88 );
1968     case('B'): return( 37 );
1969     case('E'): return( 75 );
1970     case('H'): return( 45 );
1971     case('N'): return( 86 );
1972     case('U'): return( 44 );
1973     }
1974     break;
1975    
1976     case('S'): switch( ch )
1977     {
1978     case('B'): return( 51 );
1979     case('C'): return( 21 );
1980     case('E'): return( 34 );
1981     case('I'): return( 14 );
1982     case('M'): return( 62 );
1983     case('N'): return( 50 );
1984     case('R'): return( 38 );
1985     case(' '): return( 16 );
1986     }
1987     break;
1988    
1989     case('T'): switch( ch )
1990     {
1991     case('A'): return( 73 );
1992     case('B'): return( 65 );
1993     case('C'): return( 43 );
1994     case('E'): return( 52 );
1995     case('H'): return( 90 );
1996     case('I'): return( 22 );
1997     case('L'): return( 81 );
1998     case('M'): return( 69 );
1999     }
2000     break;
2001    
2002     case('U'): if( ch==' ' )
2003     return( 92 );
2004     break;
2005    
2006     case('V'): if( ch==' ' )
2007     return( 23 );
2008     break;
2009    
2010     case('W'): if( ch==' ' )
2011     return( 74 );
2012     break;
2013    
2014     case('X'): if( ch=='E' )
2015     return( 54 );
2016     break;
2017    
2018     case('Y'): if( ch=='B' )
2019     {
2020     return( 70 );
2021     }
2022     else if( ch==' ' )
2023     return( 39 );
2024     break;
2025    
2026     case('Z'): if( ch=='N' )
2027     {
2028     return( 30 );
2029     }
2030     else if( ch=='R' )
2031     return( 40 );
2032     break;
2033     }
2034    
2035     if( (*ptr>='0') && (*ptr<='9') )
2036     if( (ch=='H') || (ch=='D') )
2037     return( 1 ); /* Hydrogen */
2038    
2039     return( 0 );
2040     }
2041    
2042     char *OBChainsParser::ParseSmiles(char *ptr, int prev)
2043     {
2044     char *name;
2045     int atomid;
2046     int next;
2047     int type;
2048     int ch;
2049    
2050     type = 0;
2051     while( (ch = *ptr++) )
2052     {
2053     switch( ch )
2054     {
2055     case('-'): type = BF_SINGLE;
2056     break;
2057     case('='): type = BF_DOUBLE;
2058     break;
2059     case('#'): type = BF_TRIPLE;
2060     break;
2061     case('^'): type = BF_SINGLE|BF_AROMATIC;
2062     break;
2063     case('~'): type = BF_DOUBLE|BF_AROMATIC;
2064     break;
2065    
2066     case(')'): return( ptr );
2067     case('.'): prev = -1;
2068     break;
2069     case('('): ptr = ParseSmiles(ptr,prev);
2070     break;
2071    
2072     default:
2073     atomid = ch-'0';
2074     while( isdigit(*ptr) )
2075     atomid = (atomid*10)+(*ptr++)-'0';
2076    
2077     for( next=0; next<MonoAtomCount; next++ )
2078     if( MonoAtom[next].atomid == atomid )
2079     break;
2080    
2081     if( next == MonoAtomCount )
2082     {
2083     name = ChainsAtomName[atomid];
2084     MonoAtom[next].elem = IdentifyElement(name);
2085     MonoAtom[next].atomid = atomid;
2086     MonoAtom[next].bcount = 0;
2087     MonoAtomCount++;
2088     }
2089    
2090     if( prev != -1 )
2091     {
2092     MonoBond[MonoBondCount].flag = type;
2093     MonoBond[MonoBondCount].src = prev;
2094     MonoBond[MonoBondCount].dst = next;
2095     MonoBondCount++;
2096    
2097     MonoAtom[prev].bcount++;
2098     MonoAtom[next].bcount++;
2099     }
2100     prev = next;
2101     }
2102     }
2103     return( ptr-1 );
2104     }
2105    
2106     #ifdef _I_WANT_TO_OUTPUT_PDB_
2107    
2108     static ChainsAtom *PDBOrder[MaxChainsAtom];
2109    
2110     int PDBSort(ChainsAtom **arg1, ChainsAtom **arg2)
2111     {
2112     ChainsAtom *atom1;
2113     ChainsAtom *atom2;
2114    
2115     atom1 = *arg1;
2116     atom2 = *arg2;
2117    
2118     if( atom1->chain != atom2->chain )
2119     return( atom1->chain - atom2->chain );
2120    
2121     if( atom1->hetflag != atom2->hetflag )
2122     return( atom1->hetflag? 1 : -1 );
2123    
2124     if( atom1->resno != atom2->resno )
2125     return( atom1->resno - atom2->resno );
2126    
2127     if( (atom1->elem==1) && (atom2->elem!=1) )
2128     return( 1 );
2129     if( (atom1->elem!=1) && (atom2->elem==1) )
2130     return( -1 );
2131    
2132     if( atom1->atomid != atom2->atomid )
2133     return( atom1->atomid - atom2->atomid );
2134    
2135     if( (atom1->elem==1) && (atom2->elem==1) )
2136     return( atom1->hcount - atom2->hcount );
2137     return( 0 );
2138     }
2139    
2140     static void OutputPDBFile(ChainsMolecule *mol, FILE *fp)
2141     {
2142     int src,dst;
2143     ChainsAtom *atom;
2144     char *ptr;
2145     int i;
2146    
2147     for( i=0; i<mol->acount; i++ )
2148     PDBOrder[i] = &mol->atom[i];
2149    
2150     #ifdef __STDC__
2151    
2152     qsort(PDBOrder,mol->acount,sizeof(ChainsAtom*),
2153     (int(*)(const void*,const void*))PDBSort);
2154     #else
2155    
2156     qsort(PDBOrder,mol->acount,sizeof(ChainsAtom*),PDBSort);
2157     #endif
2158    
2159     ptr = mol->name;
2160     while( *ptr == ' ' )
2161     ptr++;
2162    
2163     if( *ptr )
2164     {
2165     fputs("COMPND ",fp);
2166     while( *ptr )
2167     fputc(*ptr++,fp);
2168     fputc('\n',fp);
2169     }
2170    
2171     for( i=0; i<mol->acount; i++ )
2172     {
2173     atom = PDBOrder[i];
2174     atom->serno = i+1;
2175    
2176     if( atom->hetflag )
2177     {
2178     fputs("HETATM ",fp);
2179     }
2180     else
2181     fputs("ATOM ",fp);
2182    
2183     fprintf(fp,"%4d ",atom->serno);
2184    
2185     if( atom->atomid == -1 )
2186     {
2187     fprintf(fp,"%s ", etab.GetSymbol(atom->elem));
2188     }
2189     else if( atom->elem == 1 )
2190     {
2191     if( atom->hcount )
2192     {
2193     fputc(atom->hcount+'0',fp);
2194     }
2195     else
2196     fputc(' ',fp);
2197     fprintf(fp,"H%.2s",ChainsAtomName[atom->atomid]+2);
2198     }
2199     else
2200     fprintf(fp,"%.4s",ChainsAtomName[atom->atomid]);
2201    
2202     fprintf(fp," %s ",ChainsResName[atom->resid]);
2203     fprintf(fp,"%c%4d",atom->chain,atom->resno);
2204     fprintf(fp," %8.3lf%8.3lf%8.3lf",atom->x,atom->y,atom->z);
2205     fputs(" 1.00 0.00\n",fp);
2206     }
2207    
2208     for( i=0; i<mol->bcount; i++ )
2209     if( mol->bond[i].flag & BF_DOUBLE )
2210     {
2211     src = mol->atom[mol->bond[i].src].serno;
2212     dst = mol->atom[mol->bond[i].dst].serno;
2213     fprintf(fp,"CONECT%5d%5d%5d\n",src,dst,dst);
2214     fprintf(fp,"CONECT%5d%5d%5d\n",dst,src,src);
2215     }
2216     fputs("END \n",fp);
2217     }
2218    
2219     #endif
2220    
2221     } // end namespace OpenBabel
2222    
2223     //! \file chains.cpp
2224     //! \brief Parse for macromolecule chains and residues.