1import three;
2import cpkcolors;
3
4// A sample Protein Data Bank file for this example is available from
5// http://ndbserver.rutgers.edu/files/ftp/NDB/coordinates/na-biol/100d.pdb1
6
7currentlight=White;
8//currentlight=nolight;
9
10defaultrender.merge=true;  // Fast low-quality rendering
11//defaultrender.merge=false; // Slow high-quality rendering
12bool pixel=false; // Set to true to draw dots as pixels.
13real width=6;
14
15size(200);
16currentprojection=perspective(30,30,15);
17
18pen chainpen=green;
19pen hetpen=purple;
20
21string filename="100d.pdb1";
22//string filename=getstring("filename");
23
24string prefix=stripextension(filename);
25file data=input(filename);
26
27pen color(string e)
28{
29  e=replace(e," ","");
30  int n=length(e);
31  if(n < 1) return currentpen;
32  if(n > 1) e=substr(e,0,1)+downcase(substr(e,1,n-1));
33  int index=find(Element == e);
34  if(index < 0) return currentpen;
35  return rgb(Hexcolor[index]);
36}
37
38// ATOM
39string[] name,altLoc,resName,chainID,iCode,element,charge;
40int[] serial,resSeq;
41real[][] occupancy,tempFactor;
42
43bool newchain=true;
44
45struct bond
46{
47  int i,j;
48  void operator init(int i, int j) {
49    this.i=i;
50    this.j=j;
51  }
52}
53
54bond[] bonds;
55
56struct atom
57{
58  string name;
59  triple v;
60  void operator init(string name, triple v) {
61    this.name=name;
62    this.v=v;
63  }
64}
65
66struct chain
67{
68  int[] serial;
69  atom[] a;
70}
71
72int[] serials;
73chain[] chains;
74atom[] atoms;
75
76while(true) {
77  string line=data;
78  if(eof(data)) break;
79  string record=replace(substr(line,0,6)," ","");
80  if(record == "TER") {newchain=true; continue;}
81  bool ATOM=record == "ATOM";
82  bool HETATOM=record == "HETATM";
83  int serial;
84
85  atom a;
86  if(ATOM || HETATOM) {
87    serial=(int) substr(line,6,5);
88    a.name=substr(line,76,2);
89    a.v=((real) substr(line,30,8),
90         (real) substr(line,38,8),
91         (real) substr(line,46,8));
92  }
93  if(ATOM) {
94    if(newchain) {
95      chains.push(new chain);
96      newchain=false;
97    }
98    chain c=chains[chains.length-1];
99    c.serial.push(serial);
100    c.a.push(a);
101    continue;
102  }
103  if(HETATOM) {
104    serials.push(serial);
105    atoms.push(a);
106  }
107  if(record == "CONECT") {
108    int k=0;
109    int i=(int) substr(line,6,5);
110    while(true) {
111      string s=replace(substr(line,11+k,5)," ","");
112      if(s == "") break;
113      k += 5;
114      int j=(int) s;
115      if(j <= i) continue;
116      bonds.push(bond(i,j));
117    }
118  }
119}
120
121write("Number of atomic chains: ",chains.length);
122
123int natoms;
124begingroup3("chained");
125for(chain c : chains) {
126  for(int i=0; i < c.a.length-1; ++i)
127    draw(c.a[i].v--c.a[i+1].v,chainpen,currentlight);
128  for(atom a : c.a)
129    if(pixel)
130      pixel(a.v,color(a.name),width);
131    else
132      dot(a.v,color(a.name),currentlight);
133  natoms += c.a.length;
134}
135endgroup3();
136
137write("Number of chained atoms: ",natoms);
138write("Number of hetero atoms: ",atoms.length);
139
140begingroup3("hetero");
141for(atom h : atoms)
142  if(pixel)
143    pixel(h.v,color(h.name),width);
144  else
145    dot(h.v,color(h.name),currentlight);
146endgroup3();
147
148write("Number of hetero bonds: ",bonds.length);
149
150begingroup3("bonds");
151for(bond b : bonds) {
152  triple v(int i) {return atoms[find(serials == i)].v;}
153  draw(v(b.i)--v(b.j),hetpen,currentlight);
154}
155endgroup3();
156
157string options;
158string viewfilename=prefix+".views";
159
160if(!error(input(viewfilename,check=false)))
161  options="3Dviews="+locatefile(viewfilename);
162
163shipout(options=options);
164