1 /* Copyright 2004 2 Stanford University 3 4 This file is part of the DSR PDB Library. 5 6 The DSR PDB Library is free software; you can redistribute it and/or modify 7 it under the terms of the GNU Lesser General Public License as published by 8 the Free Software Foundation; either version 2.1 of the License, or (at your 9 option) any later version. 10 11 The DSR PDB Library is distributed in the hope that it will be useful, but 12 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 13 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 14 License for more details. 15 16 You should have received a copy of the GNU Lesser General Public 17 License along with the DSR PDB Library; see the file COPYING.LIB. If 18 not, write to the Free Software Foundation, Inc., 59 Temple Place - 19 Suite 330, Boston, MA 02111-1307, USA. */ 20 21 #include "dsrpdb/Residue.h" 22 #include "Residue_data.h" 23 #include <algorithm> 24 #include <iostream> 25 #include <cassert> 26 #include <dsrpdb_internal/Error_logger.h> 27 #include <limits> 28 #include <sstream> 29 30 namespace dsrpdb { 31 32 33 static Atom dummy_atom_; 34 set_has_bonds(bool tf)35 void Residue::set_has_bonds(bool tf) { 36 typedef Residue_data::Possible_bond Possible_bond; 37 38 if (!tf) { 39 bonds_.clear(); 40 } else { 41 const std::vector<Possible_bond >& bonds= Residue_data::amino_acid_data_[type()].bonds; 42 for (unsigned int i=0; i< bonds.size(); ++i){ 43 Const_atoms_iterator itf= atoms_.find(bonds[i].first); 44 Const_atoms_iterator its= atoms_.find(bonds[i].second); 45 if ( itf != atoms_.end() && its != atoms_.end()){ 46 bonds_.push_back(Bond(itf->second.index(), its->second.index())); 47 } 48 } 49 } 50 } 51 52 53 54 55 56 has_atom(Residue::Atom_label ial) const57 bool Residue::has_atom(Residue::Atom_label ial) const { 58 Residue::Atom_label al= Residue_data::fix_atom_label(label_, ial); 59 assert(can_have_atom(al)); 60 return atoms_.find(al) != atoms_.end(); 61 } 62 63 64 65 66 67 68 69 70 can_have_atom(Residue::Atom_label ial) const71 bool Residue::can_have_atom(Residue::Atom_label ial) const { 72 if (ial== AL_INVALID) return false; // hack to avoid some circularity 73 Residue::Atom_label al= Residue_data::fix_atom_label(label_, ial); 74 //assert(Residue_data::amino_acid_data_.find(label_)!= Residue_data::amino_acid_data_.end()); 75 for (unsigned int i=0; i< Residue_data::amino_acid_data_[label_].atoms.size(); ++i){ 76 if (Residue_data::amino_acid_data_[label_].atoms[i] == al) return true; 77 } 78 return false; 79 } 80 81 82 83 84 85 86 87 atom(Residue::Atom_label al) const88 const Atom &Residue::atom(Residue::Atom_label al) const { 89 Residue::Atom_label fal= Residue_data::fix_atom_label(label_, al); 90 //int fa= find_atom(al); 91 if (atoms_.find(fal) != atoms_.end()) return atoms_.find(fal)->second; 92 else { 93 return dummy_atom_; 94 } 95 } 96 97 98 99 100 101 set_atom(Residue::Atom_label ial,const Atom & a)102 void Residue::set_atom(Residue::Atom_label ial, const Atom &a) { 103 Residue::Atom_label al= Residue_data::fix_atom_label(label_, ial); 104 if (!can_have_atom(al)){ 105 dsrpdb_internal::error_logger.new_warning((std::string("Trying to set invalid atom ") + atom_label_string(ial) 106 + " on a residue of type " + type_string(label_)).c_str()); 107 } 108 if (al == AL_INVALID) { 109 return; 110 } 111 //assert(atoms_.find(al)== atoms_.end()); 112 //bonds_.clear(); 113 atoms_[al]=a; 114 // This is a bit evil, I haven't figured out a better way though. 115 atoms_[al].set_type(element(al)); 116 if (min_atom_index_) { 117 min_atom_index_= std::min(min_atom_index_, a.index()); 118 } else { 119 min_atom_index_= a.index(); 120 } 121 if (has_bonds()){ 122 set_has_bonds(false); 123 set_has_bonds(true); 124 } 125 //atoms_.push_back(a); 126 } 127 128 129 130 131 132 133 134 type() const135 Residue::Type Residue::type() const { 136 return label_; 137 } 138 dump(std::ostream & out) const139 void Residue::dump(std::ostream &out) const { 140 out << "Type: " << type_string(type()) << std::endl; 141 const std::vector<Residue::Atom_label> &valid_atoms= Residue_data::amino_acid_data_[type()].atoms; 142 for (unsigned int i=0; i< valid_atoms.size(); ++i){ 143 Residue::Atom_label al= valid_atoms[i]; 144 out << Residue::atom_label_string(al); //Residue::write_atom_label(al, out); 145 if (has_atom(al)){ 146 out << " (" << atom(al).cartesian_coords() << ") " << atom(al).index() << std::endl; 147 } else { 148 out << " X" << std::endl; 149 } 150 } 151 } 152 153 struct Compare_index{ operator ()dsrpdb::Compare_index154 bool operator()(const std::pair<Residue::Atom_label, Atom> &a, 155 const std::pair<Residue::Atom_label, Atom>&b) const { 156 return a.second.index() < b.second.index(); 157 } 158 }; 159 write(char chain,std::ostream & out) const160 void Residue::write(char chain, std::ostream &out) const { 161 static const char atom_line_oformat_[]= 162 "ATOM %5d %4s%1c%3s %1c%4d%1c %8.3f%8.3f%8.3f%6.2f%6.2f %4s%2s%2s"; 163 char line[81]; 164 std::vector<std::pair<Atom_label, Atom> > atoms(atoms_.begin(), atoms_.end()); 165 std::sort(atoms.begin(), atoms.end(), Compare_index()); 166 167 for (unsigned int i=0; i< atoms.size(); ++i) { 168 Atom_label al= atoms[i].first; 169 //Point pt= res->cartesian_coords(al); 170 const Atom &a= atoms[i].second; 171 Point pt = a.cartesian_coords(); 172 char alt=' '; 173 //char chain=' '; 174 char insertion_residue_code=' '; 175 sprintf(line, atom_line_oformat_, 176 static_cast<unsigned int>(a.index()), Residue::atom_label_string(al).c_str(), alt, 177 Residue::type_string(type()).c_str(), chain, static_cast<unsigned int>(index()), insertion_residue_code, 178 pt.x(), pt.y(), pt.z(), a.occupancy(), a.temperature_factor(), a.segment_id(), 179 a.element(), a.charge()); 180 out << line << std::endl; 181 //++anum; 182 } 183 } 184 last_atom_index() const185 Atom::Index Residue::last_atom_index() const { 186 Atom::Index max= atoms_.begin()->second.index(); 187 for ( Const_atoms_iterator it= atoms_.begin(); it != atoms_.end(); 188 ++it) { 189 max= std::max(max, it->second.index()); 190 } 191 192 return max; 193 } 194 number_of_atoms() const195 unsigned int Residue::number_of_atoms() const { 196 return atoms_.size(); 197 } 198 number_of_bonds() const199 unsigned int Residue::number_of_bonds() const { 200 bonds_begin(); 201 return bonds_.size(); 202 } 203 set_index(Residue::Index i)204 void Residue::set_index(Residue::Index i) { 205 index_=i; 206 } 207 208 atom_label(Atom::Index model_index) const209 Residue::Atom_label Residue::atom_label(Atom::Index model_index) const{ 210 Const_atoms_iterator it= atoms_iterator_from_index(model_index); 211 if (it != atoms_end()) { 212 return it->first; 213 } else { 214 return AL_INVALID; 215 } 216 } 217 218 219 220 // protected index(Residue::Atom_label al) const221 Atom::Index Residue::index(Residue::Atom_label al) const { 222 return atom(al).index(); 223 } 224 atoms_iterator_from_index(Atom::Index ind)225 Residue::Atoms_iterator Residue::atoms_iterator_from_index(Atom::Index ind) { 226 for (Atoms_iterator it= atoms_begin(); it != atoms_end(); ++it){ 227 if (it->second.index()==ind) return it; 228 } 229 dsrpdb_internal::error_logger.new_warning("Invalid atom index used to request atom from residue."); 230 231 return atoms_end(); 232 } atoms_iterator_from_index(Atom::Index ind) const233 Residue::Const_atoms_iterator Residue::atoms_iterator_from_index(Atom::Index ind) const { 234 for (Const_atoms_iterator it= atoms_begin(); it != atoms_end(); ++it){ 235 if (it->second.index()==ind) return it; 236 } 237 dsrpdb_internal::error_logger.new_warning("Invalid atom index used to request atom from residue."); 238 239 return atoms_end(); 240 } 241 242 243 244 Residue(Type al)245 Residue::Residue(Type al): atoms_(20), label_(al){ 246 Residue_data::initialize(); 247 assert(al!= INV); 248 } 249 250 sidechain_point() const251 Point Residue::sidechain_point() const { 252 double x=0; 253 double y=0; 254 double z=0; 255 int count=std::distance(Residue_data::amino_acid_data_[type()].extremes.begin(), 256 Residue_data::amino_acid_data_[type()].extremes.end()); 257 for (std::vector<Atom_label>::const_iterator it = Residue_data::amino_acid_data_[type()].extremes.begin(); 258 it != Residue_data::amino_acid_data_[type()].extremes.end(); ++it){ 259 x+= atom(*it).cartesian_coords().x(); 260 y+= atom(*it).cartesian_coords().y(); 261 z+= atom(*it).cartesian_coords().z(); 262 } 263 if (count == 0) { 264 return atom(AL_CA).cartesian_coords(); 265 } else { 266 return Point(x/count, y/count, z/count); 267 } 268 } 269 270 271 /*Residue::Simplified_atoms_iterator Residue::simplified_atoms_begin() const { 272 return Residue_data::amino_acid_data_[type()].extremes.begin(); 273 } 274 275 Residue::Simplified_atoms_iterator Residue::simplified_atoms_end() const{ 276 return Residue_data::amino_acid_data_[type()].extremes.end(); 277 }*/ 278 279 atom_label(const char * nm)280 Residue::Atom_label Residue::atom_label(const char *nm) { 281 Residue_data::initialize(); 282 char nn[5]; 283 sscanf(nm, "%4s", nn); 284 std::string s(nn); 285 for (unsigned int i=0; Residue_data::clean_atom_name_data_[i].l != AL_INVALID; ++i){ 286 if (s==Residue_data::clean_atom_name_data_[i].s){ 287 return Residue_data::clean_atom_name_data_[i].l; 288 } 289 } 290 dsrpdb_internal::error_logger.new_warning(std::string(s + " is not a known atom type.").c_str());; 291 return Residue::AL_OTHER; 292 } 293 294 295 296 297 element(Atom_label al)298 Atom::Type Residue::element(Atom_label al){ 299 Residue_data::initialize(); 300 for (unsigned int i=0; Residue_data::atom_name_data_[i].l != Residue::AL_INVALID; ++i){ 301 if (al==Residue_data::atom_name_data_[i].l){ 302 return Residue_data::atom_name_data_[i].t; 303 } 304 } 305 dsrpdb_internal::error_logger.new_internal_error("Unknown element label "); 306 return Atom::INVALID; 307 } 308 309 310 311 312 atom_label_string(Atom_label al)313 std::string Residue::atom_label_string(Atom_label al) { 314 Residue_data::initialize(); 315 for (unsigned int i=0; Residue_data::atom_name_data_[i].l != Residue::AL_INVALID; ++i){ 316 if (al==Residue_data::atom_name_data_[i].l){ 317 return Residue_data::atom_name_data_[i].s; 318 } 319 } 320 321 std::ostringstream oss; 322 oss << "Unknown atom label: " << al << " returning UNKN"; 323 dsrpdb_internal::error_logger.new_warning(oss.str().c_str()); 324 325 return "UNKN"; 326 } 327 328 329 330 331 332 type(const std::string & s)333 Residue::Type Residue::type(const std::string &s){ 334 /*switch(s[0]) { 335 case 'A': 336 switch(s[1]) { 337 case 'C': 338 return ACE; 339 case 'L': 340 return 341 default: 342 return INV; 343 }; 344 default: 345 return INV; 346 }*/ 347 if (s =="ACE") return ACE; 348 if (s =="ALA") return ALA; 349 if (s =="ARG") return ARG; 350 if (s =="ASN") return ASN; 351 if (s =="ASP") return ASP; 352 if (s =="CYS") return CYS; 353 if (s =="GLN") return GLN; 354 if (s =="GLU") return GLU; 355 if (s =="HIS") return HIS; 356 if (s =="ILE") return ILE; 357 if (s =="LEU") return LEU; 358 if (s =="LYS") return LYS; 359 if (s =="MET") return MET; 360 if (s =="NH2") return NH2; 361 if (s =="PHE") return PHE; 362 if (s =="PRO") return PRO; 363 if (s =="SER") return SER; 364 if (s =="THR") return THR; 365 if (s =="TRP") return TRP; 366 if (s =="TYR") return TYR; 367 if (s =="VAL") return VAL; 368 if (s =="GLY") return GLY; 369 else return INV; 370 } 371 type_string(Residue::Type rl)372 std::string Residue::type_string(Residue::Type rl){ 373 switch (rl) { 374 case GLY: return std::string("GLY"); 375 case ALA: return std::string("ALA"); 376 case VAL: return std::string("VAL"); 377 case LEU: return std::string("LEU"); 378 case ILE: return std::string("ILE"); 379 case SER: return std::string("SER"); 380 case THR: return std::string("THR"); 381 case CYS: return std::string("CYS"); 382 case MET: return std::string("MET"); 383 case PRO: return std::string("PRO"); 384 case ASP: return std::string("ASP"); 385 case ASN: return std::string("ASN"); 386 case GLU: return std::string("GLU"); 387 case GLN: return std::string("GLN"); 388 case LYS: return std::string("LYS"); 389 case ARG: return std::string("ARG"); 390 case HIS: return std::string("HIS"); 391 case PHE: return std::string("PHE"); 392 case TYR: return std::string("TYR"); 393 case TRP: return std::string("TRP"); 394 case ACE: return std::string("ACE"); 395 case NH2: return std::string("NH2"); 396 default: return std::string("UNK"); 397 } 398 399 } 400 401 402 403 }; 404