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