1 // $Id: mmdb_atom.cpp $ 2 // ================================================================= 3 // 4 // CCP4 Coordinate Library: support of coordinate-related 5 // functionality in protein crystallography applications. 6 // 7 // Copyright (C) Eugene Krissinel 2000-2015. 8 // 9 // This library is free software: you can redistribute it and/or 10 // modify it under the terms of the GNU Lesser General Public 11 // License version 3, modified in accordance with the provisions 12 // of the license to address the requirements of UK law. 13 // 14 // You should have received a copy of the modified GNU Lesser 15 // General Public License along with this library. If not, copies 16 // may be downloaded from http://www.ccp4.ac.uk/ccp4license.php 17 // 18 // This program is distributed in the hope that it will be useful, 19 // but WITHOUT ANY WARRANTY; without even the implied warranty of 20 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 21 // GNU Lesser General Public License for more details. 22 // 23 // ================================================================= 24 // 25 // 09.03.16 <-- Date of Last Modification. 26 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 27 // ----------------------------------------------------------------- 28 // 29 // **** Module : MMDB_Atom <implementation> 30 // ~~~~~~~~~ 31 // **** Project : MacroMolecular Data Base (MMDB) 32 // ~~~~~~~~~ 33 // **** Classes : mmdb::Atom ( atom class ) 34 // ~~~~~~~~~ mmdb::Residue ( residue class ) 35 // **** Functions: mmdb::BondAngle 36 // ~~~~~~~~~~ 37 // 38 // Copyright (C) E. Krissinel 2000-2016 39 // 40 // ================================================================= 41 // 42 43 #include <string.h> 44 #include <stdlib.h> 45 #include <math.h> 46 47 #include "mmdb_chain.h" 48 #include "mmdb_model.h" 49 #include "mmdb_root.h" 50 #include "mmdb_tables.h" 51 #include "mmdb_cifdefs.h" 52 #include "hybrid_36.h" 53 54 55 namespace mmdb { 56 57 // ================================================================ 58 59 #define ASET_CompactBinary 0x10000000 60 #define ASET_ShortTer 0x20000000 61 #define ASET_ShortHet 0x40000000 62 63 bool ignoreSegID = false; 64 bool ignoreElement = false; 65 bool ignoreCharge = false; 66 bool ignoreNonCoorPDBErrors = false; 67 bool ignoreUnmatch = false; 68 69 70 // ========================== Atom ============================= 71 Atom()72 Atom::Atom() : UDData() { 73 InitAtom(); 74 } 75 Atom(PResidue res)76 Atom::Atom ( PResidue res ) : UDData() { 77 InitAtom(); 78 if (res) 79 res->AddAtom ( this ); 80 } 81 Atom(io::RPStream Object)82 Atom::Atom ( io::RPStream Object ) : UDData(Object) { 83 InitAtom(); 84 } 85 ~Atom()86 Atom::~Atom() { 87 PPAtom A; 88 int nA; 89 FreeMemory(); 90 if (residue) { 91 A = NULL; 92 nA = 0; 93 if (residue->chain) { 94 if (residue->chain->model) { 95 A = residue->chain->model->GetAllAtoms(); 96 nA = residue->chain->model->GetNumberOfAllAtoms(); 97 } 98 } 99 residue->_ExcludeAtom ( index ); 100 if ((0<index) && (index<=nA)) A[index-1] = NULL; 101 } 102 } 103 InitAtom()104 void Atom::InitAtom() { 105 serNum = -1; // serial number 106 index = -1; // index in the file 107 name[0] = char(0); // atom name 108 label_atom_id[0] = char(0); // assigned atom name (not aligned) 109 altLoc[0] = char(0); // alternate location indicator 110 residue = NULL; // reference to residue 111 x = 0.0; // orthogonal x-coordinate in angstroms 112 y = 0.0; // orthogonal y-coordinate in angstroms 113 z = 0.0; // orthogonal z-coordinate in angstroms 114 occupancy = 0.0; // occupancy 115 tempFactor = 0.0; // temperature factor 116 segID[0] = char(0); // segment identifier 117 strcpy ( element," " ); // chemical element symbol - RIGHT JUSTIFIED 118 energyType[0] = char(0); // chemical element symbol - RIGHT JUSTIFIED 119 charge = 0.0; // charge on the atom 120 sigX = 0.0; // standard deviation of the stored x-coord 121 sigY = 0.0; // standard deviation of the stored y-coord 122 sigZ = 0.0; // standard deviation of the stored z-coord 123 sigOcc = 0.0; // standard deviation of occupancy 124 sigTemp = 0.0; // standard deviation of temperature factor 125 u11 = 0.0; // 126 u22 = 0.0; // anisotropic 127 u33 = 0.0; // 128 u12 = 0.0; // temperature 129 u13 = 0.0; // 130 u23 = 0.0; // factors 131 su11 = 0.0; // 132 su22 = 0.0; // standard 133 su33 = 0.0; // deviations of 134 su12 = 0.0; // anisotropic 135 su13 = 0.0; // temperature 136 su23 = 0.0; // factors 137 Het = false; // indicator of atom in non-standard groups 138 Ter = false; // chain terminator 139 WhatIsSet = 0x00000000; // nothing is set 140 nBonds = 0; // no bonds 141 Bond = NULL; // empty array of bonds 142 } 143 FreeMemory()144 void Atom::FreeMemory() { 145 FreeBonds(); 146 } 147 FreeBonds()148 void Atom::FreeBonds() { 149 if (Bond) delete[] Bond; 150 Bond = NULL; 151 nBonds = 0; 152 } 153 GetNBonds()154 int Atom::GetNBonds() { 155 return nBonds & 0x000000FF; 156 } 157 GetBonds(RPAtomBond atomBond,int & nAtomBonds)158 void Atom::GetBonds ( RPAtomBond atomBond, int & nAtomBonds ) { 159 // This GetBonds(..) returns pointer to the Atom's 160 // internal Bond structure, IT MUST NOT BE DISPOSED. 161 nAtomBonds = nBonds & 0x000000FF; 162 atomBond = Bond; 163 } 164 GetBonds(RPAtomBondI atomBondI,int & nAtomBonds)165 void Atom::GetBonds ( RPAtomBondI atomBondI, int & nAtomBonds ) { 166 // This GetBonds(..) disposes AtomBondI, if it was not set 167 // to NULL, allocates AtomBondI[nAtomBonds] and returns its 168 // pointer. AtomBondI MUST BE DISPOSED BY APPLICATION. 169 int i; 170 171 if (atomBondI) delete[] atomBondI; 172 173 nAtomBonds = nBonds & 0x000000FF; 174 175 if (nAtomBonds<=0) 176 atomBondI = NULL; 177 else { 178 atomBondI = new AtomBondI[nAtomBonds]; 179 for (i=0;i<nAtomBonds;i++) { 180 if (Bond[i].atom) 181 atomBondI[i].index = Bond[i].atom->index; 182 else atomBondI[i].index = -1; 183 atomBondI[i].order = Bond[i].order; 184 } 185 186 } 187 188 } 189 GetBonds(PAtomBondI AtomBondI,int & nAtomBonds,int maxlength)190 void Atom::GetBonds ( PAtomBondI AtomBondI, int & nAtomBonds, 191 int maxlength ) { 192 // This GetBonds(..) does not dispose or allocate AtomBond. 193 // It is assumed that length of AtomBond is sufficient to 194 // accomodate all bonded atoms. 195 int i; 196 197 nAtomBonds = IMin(maxlength,nBonds & 0x000000FF); 198 199 for (i=0;i<nAtomBonds;i++) { 200 if (Bond[i].atom) 201 AtomBondI[i].index = Bond[i].atom->index; 202 else AtomBondI[i].index = -1; 203 AtomBondI[i].order = Bond[i].order; 204 } 205 206 } 207 208 AddBond(PAtom bond_atom,int bond_order,int nAdd_bonds)209 int Atom::AddBond ( PAtom bond_atom, int bond_order, 210 int nAdd_bonds ) { 211 PAtomBond B1; 212 int i,k,nb,nballoc; 213 214 nb = nBonds & 0x000000FF; 215 k = -1; 216 for (i=0;(i<nb) && (k<0);i++) 217 if (Bond[i].atom==bond_atom) k = i; 218 if (k>=0) return -k; 219 220 nballoc = (nBonds >> 8) & 0x000000FF; 221 if (nBonds>=nballoc) { 222 nballoc += nAdd_bonds; 223 B1 = new AtomBond[nballoc]; 224 for (i=0;i<nb;i++) { 225 B1[i].atom = Bond[i].atom; 226 B1[i].order = Bond[i].order; 227 } 228 if (Bond) delete[] Bond; 229 Bond = B1; 230 } 231 Bond[nb].atom = bond_atom; 232 Bond[nb].order = bond_order; 233 nb++; 234 235 nBonds = nb | (nballoc << 8); 236 237 return nb; 238 239 } 240 SetResidue(PResidue res)241 void Atom::SetResidue ( PResidue res ) { 242 residue = res; 243 } 244 StandardPDBOut(cpstr Record,pstr S)245 void Atom::StandardPDBOut ( cpstr Record, pstr S ) { 246 char N[10]; 247 strcpy ( S,Record ); 248 PadSpaces ( S,80 ); 249 if (serNum>99999) { 250 hy36encode ( 5,serNum,N ); 251 strcpy_n ( &(S[6]),N,5 ); 252 } else if (serNum>0) 253 PutInteger ( &(S[6]),serNum,5 ); 254 else if (index<=99999) 255 PutInteger ( &(S[6]),index,5 ); 256 else { 257 hy36encode ( 5,index,N ); 258 strcpy_n ( &(S[6]),N,5 ); 259 } 260 if (!Ter) { 261 if (altLoc[0]) S[16] = altLoc[0]; 262 strcpy_n ( &(S[12]),name ,4 ); 263 strcpy_n ( &(S[72]),segID ,4 ); 264 strcpy_nr ( &(S[76]),element,2 ); 265 if (WhatIsSet & ASET_Charge) { 266 if (charge>0) sprintf ( N,"%1i+",mround(charge) ); 267 else if (charge<0) sprintf ( N,"%1i-",mround(-charge) ); 268 else strcpy ( N," " ); 269 strcpy_n ( &(S[78]),N,2 ); 270 } else 271 strcpy_n ( &(S[78])," ",2 ); 272 } 273 strcpy_nr ( &(S[17]),residue->name,3 ); 274 strcpy_nr ( &(S[20]),residue->chain->chainID,2 ); 275 if (residue->seqNum>MinInt4) { 276 if ((-999<=residue->seqNum) && (residue->seqNum<=9999)) 277 PutIntIns ( &(S[22]),residue->seqNum,4,residue->insCode ); 278 else { 279 hy36encode ( 4,residue->seqNum,N ); 280 strcpy_n ( &(S[22]),N,4 ); 281 } 282 } 283 } 284 PDBASCIIDump(io::RFile f)285 void Atom::PDBASCIIDump ( io::RFile f ) { 286 // makes the ASCII PDB ATOM, HETATM, SIGATOM, ANISOU 287 // SIGUIJ and TER lines from the class' data 288 char S[100]; 289 if (Ter) { 290 if (WhatIsSet & ASET_Coordinates) { 291 StandardPDBOut ( pstr("TER"),S ); 292 f.WriteLine ( S ); 293 } 294 } else { 295 if (WhatIsSet & ASET_Coordinates) { 296 if (Het) StandardPDBOut ( pstr("HETATM"),S ); 297 else StandardPDBOut ( pstr("ATOM") ,S ); 298 PutRealF ( &(S[30]),x,8,3 ); 299 PutRealF ( &(S[38]),y,8,3 ); 300 PutRealF ( &(S[46]),z,8,3 ); 301 if (WhatIsSet & ASET_Occupancy) 302 PutRealF ( &(S[54]),occupancy ,6,2 ); 303 if (WhatIsSet & ASET_tempFactor) 304 PutRealF ( &(S[60]),tempFactor,6,2 ); 305 f.WriteLine ( S ); 306 } 307 if (WhatIsSet & ASET_CoordSigma) { 308 StandardPDBOut ( pstr("SIGATM"),S ); 309 PutRealF ( &(S[30]),sigX,8,3 ); 310 PutRealF ( &(S[38]),sigY,8,3 ); 311 PutRealF ( &(S[46]),sigZ,8,3 ); 312 if ((WhatIsSet & ASET_OccSigma) && 313 (WhatIsSet & ASET_Occupancy)) 314 PutRealF ( &(S[54]),sigOcc,6,2 ); 315 if ((WhatIsSet & ASET_tFacSigma) && 316 (WhatIsSet & ASET_tempFactor)) 317 PutRealF ( &(S[60]),sigTemp,6,2 ); 318 f.WriteLine ( S ); 319 } 320 if (WhatIsSet & ASET_Anis_tFac) { 321 StandardPDBOut ( pstr("ANISOU"),S ); 322 PutInteger ( &(S[28]),mround(u11*1.0e4),7 ); 323 PutInteger ( &(S[35]),mround(u22*1.0e4),7 ); 324 PutInteger ( &(S[42]),mround(u33*1.0e4),7 ); 325 PutInteger ( &(S[49]),mround(u12*1.0e4),7 ); 326 PutInteger ( &(S[56]),mround(u13*1.0e4),7 ); 327 PutInteger ( &(S[63]),mround(u23*1.0e4),7 ); 328 f.WriteLine ( S ); 329 if (WhatIsSet & ASET_Anis_tFSigma) { 330 StandardPDBOut ( pstr("SIGUIJ"),S ); 331 PutInteger ( &(S[28]),mround(su11*1.0e4),7 ); 332 PutInteger ( &(S[35]),mround(su22*1.0e4),7 ); 333 PutInteger ( &(S[42]),mround(su33*1.0e4),7 ); 334 PutInteger ( &(S[49]),mround(su12*1.0e4),7 ); 335 PutInteger ( &(S[56]),mround(su13*1.0e4),7 ); 336 PutInteger ( &(S[63]),mround(su23*1.0e4),7 ); 337 f.WriteLine ( S ); 338 } 339 } 340 } 341 } 342 MakeCIF(mmcif::PData CIF)343 void Atom::MakeCIF ( mmcif::PData CIF ) { 344 mmcif::PLoop Loop; 345 AtomName AtName; 346 Element el; 347 char N[10]; 348 int i,j,RC; 349 PChain chain = NULL; 350 PModel model = NULL; 351 //bool singleModel = true; 352 353 if (residue) chain = residue->chain; 354 if (chain) model = PModel(chain->model); 355 // if (model) { 356 // if (model->manager) 357 // singleModel = PRoot(model->manager)->nModels<=1; 358 // } 359 360 /* 361 362 loop_ 363 *0 _atom_site.group_PDB 364 *1 _atom_site.id 365 *2 _atom_site.type_symbol <- chem elem 366 -3 _atom_site.label_atom_id <- atom name 367 *4 _atom_site.label_alt_id <- alt code 368 =5 _atom_site.label_comp_id <- res name ??? 369 =6 _atom_site.label_asym_id <- chain id ??? 370 =7 _atom_site.label_entity_id < ??? 371 =8 _atom_site.label_seq_id <- poly seq id 372 +9 _atom_site.pdbx_PDB_ins_code <- ins code 373 -10 _atom_site.segment_id <- segment id 374 *11 _atom_site.Cartn_x 375 *12 _atom_site.Cartn_y 376 *13 _atom_site.Cartn_z 377 *14 _atom_site.occupancy 378 *15 _atom_site.B_iso_or_equiv 379 *16 _atom_site.Cartn_x_esd 380 *17 _atom_site.Cartn_y_esd 381 *18 _atom_site.Cartn_z_esd 382 *19 _atom_site.occupancy_esd 383 *20 _atom_site.B_iso_or_equiv_esd 384 *21 _atom_site.pdbx_formal_charge 385 +22 _atom_site.auth_seq_id <- seq id we want 386 +23 _atom_site.auth_comp_id <- res name we want 387 +24 _atom_site.auth_asym_id <- ch id we want ? 388 *25 _atom_site.auth_atom_id <- atom name we want ? 389 +26 _atom_site.pdbx_PDB_model_num <- model no 390 391 '+' is read in Root::CheckAtomPlace() 392 '=' new in residue 393 '-' new in atom 394 395 396 */ 397 398 RC = CIF->AddLoop ( CIFCAT_ATOM_SITE,Loop ); 399 if (RC!=mmcif::CIFRC_Ok) { 400 // the category was (re)created, provide tags 401 402 Loop->AddLoopTag ( CIFTAG_GROUP_PDB ); // ATOM, TER etc. 403 Loop->AddLoopTag ( CIFTAG_ID ); // serial number 404 405 Loop->AddLoopTag ( CIFTAG_TYPE_SYMBOL ); // element symbol 406 Loop->AddLoopTag ( CIFTAG_LABEL_ATOM_ID ); // atom name 407 Loop->AddLoopTag ( CIFTAG_LABEL_ALT_ID ); // alt location 408 Loop->AddLoopTag ( CIFTAG_LABEL_COMP_ID ); // residue name 409 Loop->AddLoopTag ( CIFTAG_LABEL_ASYM_ID ); // chain ID 410 Loop->AddLoopTag ( CIFTAG_LABEL_ENTITY_ID ); // entity ID 411 Loop->AddLoopTag ( CIFTAG_LABEL_SEQ_ID ); // res seq number 412 Loop->AddLoopTag ( CIFTAG_PDBX_PDB_INS_CODE ); // insertion code 413 Loop->AddLoopTag ( CIFTAG_SEGMENT_ID ); // segment ID 414 415 Loop->AddLoopTag ( CIFTAG_CARTN_X ); // x-coordinate 416 Loop->AddLoopTag ( CIFTAG_CARTN_Y ); // y-coordinate 417 Loop->AddLoopTag ( CIFTAG_CARTN_Z ); // z-coordinate 418 Loop->AddLoopTag ( CIFTAG_OCCUPANCY ); // occupancy 419 Loop->AddLoopTag ( CIFTAG_B_ISO_OR_EQUIV ); // temp factor 420 421 Loop->AddLoopTag ( CIFTAG_CARTN_X_ESD ); // x-sigma 422 Loop->AddLoopTag ( CIFTAG_CARTN_Y_ESD ); // y-sigma 423 Loop->AddLoopTag ( CIFTAG_CARTN_Z_ESD ); // z-sigma 424 Loop->AddLoopTag ( CIFTAG_OCCUPANCY_ESD ); // occupancy-sigma 425 Loop->AddLoopTag ( CIFTAG_B_ISO_OR_EQUIV_ESD ); // t-factor-sigma 426 427 Loop->AddLoopTag ( CIFTAG_PDBX_FORMAL_CHARGE ); // charge on atom 428 429 Loop->AddLoopTag ( CIFTAG_AUTH_SEQ_ID ); // res seq number 430 Loop->AddLoopTag ( CIFTAG_AUTH_COMP_ID ); // residue name 431 Loop->AddLoopTag ( CIFTAG_AUTH_ASYM_ID ); // chain id 432 Loop->AddLoopTag ( CIFTAG_AUTH_ATOM_ID ); // atom name 433 434 Loop->AddLoopTag ( CIFTAG_PDBX_PDB_MODEL_NUM ); // model number 435 436 } 437 438 /* 439 if (Ter) { // ter record 440 441 if (!(WhatIsSet & ASET_Coordinates)) 442 return; 443 444 // (0) 445 Loop->AddString ( pstr("TER") ); 446 // (1) 447 if (serNum>0) Loop->AddInteger ( serNum ); 448 else Loop->AddInteger ( index ); 449 // (2,3,4) 450 Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); // no element symbol 451 Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); // no atom name 452 Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); // no alt code 453 if (residue) { 454 // (5) 455 Loop->AddString ( residue->label_comp_id ); 456 // (6) 457 Loop->AddString ( residue->label_asym_id ); 458 // (7) 459 if (residue->label_entity_id>0) 460 Loop->AddInteger ( residue->label_entity_id ); 461 else Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 462 // (8) 463 if (residue->label_seq_id>MinInt4) 464 Loop->AddInteger ( residue->label_seq_id ); 465 else Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 466 // (9) 467 Loop->AddString ( residue->insCode,true ); 468 } else { 469 // (5,6,7,8,9) 470 Loop->AddString ( NULL ); 471 Loop->AddString ( NULL ); 472 Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 473 Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 474 Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 475 } 476 477 // (10-21) 478 for (i=10;i<=21;i++) 479 Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); 480 481 // (22,23) 482 if (residue) { 483 if (residue->seqNum>MinInt4) 484 Loop->AddInteger ( residue->seqNum ); 485 else Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 486 Loop->AddString ( residue->name ); 487 } else { 488 Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 489 Loop->AddString ( NULL ); 490 } 491 492 // (24) 493 if (chain) Loop->AddString ( chain->chainID,true ); 494 else Loop->AddString ( NULL ); 495 496 // (25) 497 Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); // no atom name 498 499 } else 500 */ 501 502 if (WhatIsSet & (ASET_Coordinates | ASET_CoordSigma)) { 503 // normal atom record 504 505 if (!WhatIsSet) return; 506 507 // (0) 508 if (Het) Loop->AddString ( pstr("HETATM") ); 509 else Loop->AddString ( pstr("ATOM") ); 510 // (1) 511 if (serNum>0) Loop->AddInteger ( serNum ); 512 else Loop->AddInteger ( index ); 513 514 515 if (WhatIsSet & ASET_Coordinates) { 516 517 // (2) 518 strcpy_css ( el,element ); 519 Loop->AddString ( el,true ); 520 // (3) 521 strcpy_css ( AtName,label_atom_id ); 522 Loop->AddString ( AtName ); // assigned atom name 523 // (4) 524 Loop->AddString ( altLoc,true ); // alt code 525 526 if (residue) { 527 // (5) 528 Loop->AddString ( residue->label_comp_id ); 529 // (6) 530 Loop->AddString ( residue->label_asym_id ); 531 // (7) 532 if (residue->label_entity_id>0) 533 Loop->AddInteger ( residue->label_entity_id ); 534 else Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 535 // (8) 536 if (residue->label_seq_id>MinInt4) 537 Loop->AddInteger ( residue->label_seq_id ); 538 else Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 539 // (9) 540 Loop->AddString ( residue->insCode,true ); 541 } else { 542 // (5,6,7,8,9) 543 Loop->AddString ( NULL ); 544 Loop->AddString ( NULL ); 545 Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 546 Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 547 Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 548 } 549 550 // (10) 551 Loop->AddString ( segID ,true ); 552 // (11,12,13) 553 Loop->AddReal ( x ); 554 Loop->AddReal ( y ); 555 Loop->AddReal ( z ); 556 // (14) 557 if (WhatIsSet & ASET_Occupancy) 558 Loop->AddReal ( occupancy ); 559 else Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); 560 // (15) 561 if (WhatIsSet & ASET_tempFactor) 562 Loop->AddReal ( tempFactor ); 563 else Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); 564 565 // (16,17,18) 566 if (WhatIsSet & ASET_CoordSigma) { 567 Loop->AddReal ( sigX ); 568 Loop->AddReal ( sigY ); 569 Loop->AddReal ( sigZ ); 570 } else { 571 Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); 572 Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); 573 Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); 574 } 575 // (19) 576 if ((WhatIsSet & ASET_OccSigma) && (WhatIsSet & ASET_Occupancy)) 577 Loop->AddReal ( sigOcc ); 578 else Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); 579 // (20) 580 if ((WhatIsSet & ASET_tFacSigma) && (WhatIsSet & ASET_tempFactor)) 581 Loop->AddReal ( sigTemp ); 582 else Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); 583 584 } else 585 for (i=0;i<18;i++) 586 Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); 587 588 // (21) 589 if (WhatIsSet & ASET_Charge) { 590 sprintf ( N,"%+2i",mround(charge) ); 591 Loop->AddString ( N,true ); 592 } else 593 Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); 594 595 if (residue) { 596 // (22) 597 if (residue->seqNum>MinInt4) 598 Loop->AddInteger ( residue->seqNum ); 599 else Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 600 // (23) 601 Loop->AddString ( residue->name ); 602 } else { 603 Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 604 Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 605 } 606 607 // (24) 608 if (chain) Loop->AddString ( chain->chainID,true ); 609 else Loop->AddNoData ( mmcif::CIF_NODATA_DOT ); 610 // (25) 611 strcpy_css ( AtName,name ); 612 Loop->AddString ( AtName ); // atom name 613 614 } 615 616 // (26) 617 if (!model) Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); 618 else if (model->serNum>0) Loop->AddInteger ( model->serNum ); 619 else Loop->AddNoData ( mmcif::CIF_NODATA_QUESTION ); 620 621 622 if (WhatIsSet & ASET_Anis_tFac) { 623 624 RC = CIF->AddLoop ( CIFCAT_ATOM_SITE_ANISOTROP,Loop ); 625 if (RC!=mmcif::CIFRC_Ok) { 626 // the category was (re)created, provide tags 627 Loop->AddLoopTag ( CIFTAG_ID ); // serial number 628 Loop->AddLoopTag ( CIFTAG_U11 ); // component u11 629 Loop->AddLoopTag ( CIFTAG_U22 ); // component u22 630 Loop->AddLoopTag ( CIFTAG_U33 ); // component u33 631 Loop->AddLoopTag ( CIFTAG_U12 ); // component u12 632 Loop->AddLoopTag ( CIFTAG_U13 ); // component u13 633 Loop->AddLoopTag ( CIFTAG_U23 ); // component u23 634 Loop->AddLoopTag ( CIFTAG_U11_ESD ); // component u11 sigma 635 Loop->AddLoopTag ( CIFTAG_U22_ESD ); // component u22 sigma 636 Loop->AddLoopTag ( CIFTAG_U33_ESD ); // component u33 sigma 637 Loop->AddLoopTag ( CIFTAG_U12_ESD ); // component u12 sigma 638 Loop->AddLoopTag ( CIFTAG_U13_ESD ); // component u13 sigma 639 Loop->AddLoopTag ( CIFTAG_U23_ESD ); // component u23 sigma 640 for (i=1;i<index;i++) { 641 Loop->AddInteger ( i ); 642 for (j=0;j<12;j++) 643 Loop->AddString ( NULL ); 644 } 645 } 646 647 if (serNum>0) Loop->AddInteger ( serNum ); 648 else Loop->AddInteger ( index ); 649 650 Loop->AddReal ( u11 ); 651 Loop->AddReal ( u22 ); 652 Loop->AddReal ( u33 ); 653 Loop->AddReal ( u12 ); 654 Loop->AddReal ( u13 ); 655 Loop->AddReal ( u23 ); 656 if (WhatIsSet & ASET_Anis_tFSigma) { 657 Loop->AddReal ( su11 ); 658 Loop->AddReal ( su22 ); 659 Loop->AddReal ( su33 ); 660 Loop->AddReal ( su12 ); 661 Loop->AddReal ( su13 ); 662 Loop->AddReal ( su23 ); 663 } 664 665 } 666 667 } 668 ConvertPDBATOM(int ix,cpstr S)669 ERROR_CODE Atom::ConvertPDBATOM ( int ix, cpstr S ) { 670 // Gets data from the PDB ASCII ATOM record. 671 // This function DOES NOT check the "ATOM" keyword and 672 // does not decode the chain and residue parameters! These 673 // must be treated by calling process, see 674 // Chain::ConvertPDBASCII(). 675 676 index = ix; 677 678 if (WhatIsSet & ASET_Coordinates) 679 return Error_ATOM_AlreadySet; 680 681 if (!(GetReal(x,&(S[30]),8) && 682 GetReal(y,&(S[38]),8) && 683 GetReal(z,&(S[46]),8))) 684 return Error_ATOM_Unrecognized; 685 686 WhatIsSet |= ASET_Coordinates; 687 Het = false; 688 Ter = false; 689 690 if (GetReal(occupancy ,&(S[54]),6)) WhatIsSet |= ASET_Occupancy; 691 if (GetReal(tempFactor,&(S[60]),6)) WhatIsSet |= ASET_tempFactor; 692 693 if (WhatIsSet & (ASET_CoordSigma | ASET_Anis_tFac | 694 ASET_Anis_tFSigma)) 695 // something was already submitted. check complience 696 return CheckData ( S ); 697 else 698 // first data submission. just take the data 699 GetData ( S ); 700 701 return Error_NoError; 702 703 } 704 SetAtomName(int ix,int sN,const AtomName aName,const AltLoc aLoc,const SegID sID,const Element eName)705 void Atom::SetAtomName ( int ix, 706 int sN, 707 const AtomName aName, 708 const AltLoc aLoc, 709 const SegID sID, 710 const Element eName ) { 711 index = ix; 712 serNum = sN; 713 strcpy ( name ,aName ); 714 strcpy ( label_atom_id,aName ); 715 strcpy_css ( altLoc ,pstr(aLoc) ); 716 strcpy_css ( segID ,pstr(sID) ); 717 if (!eName[0]) element[0] = char(0); 718 else if (!eName[1]) { 719 element[0] = ' '; 720 strcpy ( &(element[1]),eName ); 721 } else 722 strcpy ( element,eName ); 723 WhatIsSet = 0; 724 } 725 ConvertPDBSIGATM(int ix,cpstr S)726 ERROR_CODE Atom::ConvertPDBSIGATM ( int ix, cpstr S ) { 727 // Gets data from the PDB ASCII SIGATM record. 728 // This function DOES NOT check the "SIGATM" keyword and 729 // does not decode the chain and residue parameters! These 730 // must be treated by the calling process, see 731 // Chain::ConvertPDBASCII(). 732 733 index = ix; 734 735 if (WhatIsSet & ASET_CoordSigma) 736 return Error_ATOM_AlreadySet; 737 738 if (!(GetReal(sigX,&(S[30]),8) && 739 GetReal(sigY,&(S[38]),8) && 740 GetReal(sigZ,&(S[46]),8))) 741 return Error_ATOM_Unrecognized; 742 743 WhatIsSet |= ASET_CoordSigma; 744 745 if (GetReal(sigOcc ,&(S[54]),6)) WhatIsSet |= ASET_OccSigma; 746 if (GetReal(sigTemp,&(S[60]),6)) WhatIsSet |= ASET_tFacSigma; 747 748 if (WhatIsSet & (ASET_Coordinates | ASET_Anis_tFac | 749 ASET_Anis_tFSigma)) 750 // something was already submitted. check complience 751 return CheckData ( S ); 752 else 753 // first data submission. just take the data 754 GetData ( S ); 755 756 return Error_NoError; 757 758 } 759 ConvertPDBANISOU(int ix,cpstr S)760 ERROR_CODE Atom::ConvertPDBANISOU ( int ix, cpstr S ) { 761 // Gets data from the PDB ASCII ANISOU record. 762 // This function DOES NOT check the "ANISOU" keyword and 763 // does not decode chain and residue parameters! These must 764 // be treated by the calling process, see 765 // Chain::ConvertPDBASCII(). 766 767 index = ix; 768 769 if (WhatIsSet & ASET_Anis_tFac) 770 return Error_ATOM_AlreadySet; 771 772 if (!(GetReal(u11,&(S[28]),7) && 773 GetReal(u22,&(S[35]),7) && 774 GetReal(u33,&(S[42]),7) && 775 GetReal(u12,&(S[49]),7) && 776 GetReal(u13,&(S[56]),7) && 777 GetReal(u23,&(S[63]),7))) 778 return Error_ATOM_Unrecognized; 779 780 u11 /= 1.0e4; 781 u22 /= 1.0e4; 782 u33 /= 1.0e4; 783 u12 /= 1.0e4; 784 u13 /= 1.0e4; 785 u23 /= 1.0e4; 786 787 WhatIsSet |= ASET_Anis_tFac; 788 789 if (WhatIsSet & (ASET_Coordinates | ASET_CoordSigma | 790 ASET_Anis_tFSigma)) 791 // something was already submitted. check complience 792 return CheckData ( S ); 793 else 794 // first data submission. just take the data 795 GetData ( S ); 796 797 return Error_NoError; 798 799 } 800 ConvertPDBSIGUIJ(int ix,cpstr S)801 ERROR_CODE Atom::ConvertPDBSIGUIJ ( int ix, cpstr S ) { 802 // Gets data from the PDB ASCII SIGUIJ record. 803 // This function DOES NOT check the "SIGUIJ" keyword and 804 // does not decode the chain and residue parameters! These 805 // must be treated by the calling process, see 806 // Chain::ConvertPDBASCII(). 807 808 index = ix; 809 810 if (WhatIsSet & ASET_Anis_tFSigma) 811 return Error_ATOM_AlreadySet; 812 813 if (!(GetReal(su11,&(S[28]),7) && 814 GetReal(su22,&(S[35]),7) && 815 GetReal(su33,&(S[42]),7) && 816 GetReal(su12,&(S[49]),7) && 817 GetReal(su13,&(S[56]),7) && 818 GetReal(su23,&(S[63]),7))) 819 return Error_ATOM_Unrecognized; 820 821 su11 /= 1.0e4; 822 su22 /= 1.0e4; 823 su33 /= 1.0e4; 824 su12 /= 1.0e4; 825 su13 /= 1.0e4; 826 su23 /= 1.0e4; 827 828 WhatIsSet |= ASET_Anis_tFSigma; 829 830 if (WhatIsSet & (ASET_Coordinates | ASET_CoordSigma | 831 ASET_Anis_tFac)) 832 // something was already submitted. check complience 833 return CheckData ( S ); 834 else 835 // first data submission. just take the data 836 GetData ( S ); 837 838 return Error_NoError; 839 840 } 841 ConvertPDBTER(int ix,cpstr S)842 ERROR_CODE Atom::ConvertPDBTER ( int ix, cpstr S ) { 843 // Gets data from the PDB ASCII TER record. 844 // This function DOES NOT check the "TER" keyword and 845 // does not decode the chain and residue parameters! These 846 // must be treated by the calling process, see 847 // Chain::ConvertPDBASCII(). 848 849 index = ix; 850 851 if (((S[6]>='0') && (S[6]<='9')) || (S[6]==' ')) { 852 // Although against strict PDB format, 'TER' is 853 // actually allowed not to have a serial number. 854 // This negative value implies that the number is 855 // not set. 856 if (!(GetInteger(serNum,&(S[6]),5))) serNum = -1; 857 } else 858 hy36decode ( 5,&(S[6]),5,&serNum ); 859 860 // if (!(GetInteger(serNum,&(S[6]),5))) serNum = -1; 861 862 if (WhatIsSet & ASET_Coordinates) 863 return Error_ATOM_AlreadySet; 864 865 WhatIsSet |= ASET_Coordinates; 866 Het = false; 867 Ter = true; 868 name[0] = char(0); 869 label_atom_id[0] = char(0); 870 element[0] = char(0); 871 872 return Error_NoError; 873 874 } 875 876 GetModelNum()877 int Atom::GetModelNum() { 878 if (residue) { 879 if (residue->chain) 880 if (residue->chain->model) 881 return residue->chain->model->GetSerNum(); 882 } 883 return 0; 884 } 885 GetChainID()886 pstr Atom::GetChainID() { 887 if (residue) { 888 if (residue->chain) return residue->chain->chainID; 889 } 890 return pstr(""); 891 } 892 GetLabelAsymID()893 pstr Atom::GetLabelAsymID() { 894 if (residue) return residue->label_asym_id; 895 else return pstr(""); 896 } 897 GetResName()898 pstr Atom::GetResName() { 899 if (residue) return residue->name; 900 else return pstr(""); 901 } 902 GetLabelCompID()903 pstr Atom::GetLabelCompID() { 904 if (residue) return residue->label_comp_id; 905 else return pstr(""); 906 } 907 GetAASimilarity(const ResName resName)908 int Atom::GetAASimilarity ( const ResName resName ) { 909 if (residue) 910 return mmdb::GetAASimilarity ( pstr(residue->name), 911 pstr(resName) ); 912 else return -3; 913 } 914 GetAASimilarity(PAtom A)915 int Atom::GetAASimilarity ( PAtom A ) { 916 if (residue) { 917 if (A->residue) 918 return mmdb::GetAASimilarity ( pstr(residue->name), 919 pstr(A->residue->name) ); 920 else return -4; 921 } else 922 return -3; 923 } 924 GetAAHydropathy()925 realtype Atom::GetAAHydropathy() { 926 if (residue) 927 return mmdb::GetAAHydropathy ( pstr(residue->name) ); 928 else return MaxReal; 929 } 930 GetOccupancy()931 realtype Atom::GetOccupancy() { 932 if (WhatIsSet & ASET_Occupancy) return occupancy; 933 else return 0.0; 934 } 935 GetSeqNum()936 int Atom::GetSeqNum() { 937 if (residue) return residue->seqNum; 938 else return ATOM_NoSeqNum; 939 } 940 GetLabelSeqID()941 int Atom::GetLabelSeqID() { 942 if (residue) return residue->label_seq_id; 943 else return ATOM_NoSeqNum; 944 } 945 GetLabelEntityID()946 int Atom::GetLabelEntityID() { 947 if (residue) return residue->label_entity_id; 948 else return -1; 949 } 950 GetInsCode()951 pstr Atom::GetInsCode() { 952 if (residue) return residue->insCode; 953 else return pstr(""); 954 } 955 GetSSEType()956 int Atom::GetSSEType() { 957 // works only after SSE calculations 958 if (residue) return residue->SSE; 959 else return SSE_None; 960 } 961 GetAtomCharge(pstr chrg)962 pstr Atom::GetAtomCharge ( pstr chrg ) { 963 if (WhatIsSet & ASET_Charge) sprintf ( chrg,"%+2i",mround(charge) ); 964 else strcpy ( chrg," " ); 965 return chrg; 966 } 967 GetResidue()968 PResidue Atom::GetResidue() { 969 return residue; 970 } 971 GetChain()972 PChain Atom::GetChain() { 973 if (residue) return residue->chain; 974 else return NULL; 975 } 976 GetModel()977 PModel Atom::GetModel() { 978 if (residue) { 979 if (residue->chain) return (PModel)residue->chain->model; 980 } 981 return NULL; 982 } 983 GetResidueNo()984 int Atom::GetResidueNo() { 985 if (residue) { 986 if (residue->chain) 987 return residue->chain->GetResidueNo ( 988 residue->seqNum,residue->insCode ); 989 else return -2; 990 } else 991 return -1; 992 } 993 994 GetCoordHierarchy()995 void * Atom::GetCoordHierarchy() { 996 if (residue) return residue->GetCoordHierarchy(); 997 return NULL; 998 } 999 1000 GetStat(realtype v,realtype & v_min,realtype & v_max,realtype & v_m,realtype & v_m2)1001 void Atom::GetStat ( realtype v, 1002 realtype & v_min, realtype & v_max, 1003 realtype & v_m, realtype & v_m2 ) { 1004 if (v<v_min) v_min = v; 1005 if (v>v_max) v_max = v; 1006 v_m += v; 1007 v_m2 += v*v; 1008 } 1009 1010 1011 GetChainCalphas(PPAtom & Calphas,int & nCalphas,cpstr altLoc)1012 void Atom::GetChainCalphas ( PPAtom & Calphas, int & nCalphas, 1013 cpstr altLoc ) { 1014 // GetChainCalphas(...) is a specialized function for quick 1015 // access to C-alphas of chain which includes given atom. 1016 // This function works faster than an equivalent implementation 1017 // through MMDB's selection procedures. 1018 // Parameters: 1019 // Calphas - array to accept pointers on C-alpha atoms 1020 // If Calphas!=NULL, then the function will 1021 // delete and re-allocate it. When the array 1022 // is no longer needed, the application MUST 1023 // delete it: delete[] Calphas; Deleting 1024 // Calphas does not delete atoms from MMDB. 1025 // nCalphas - integer to accept number of C-alpha atoms 1026 // and the length of Calphas array. 1027 // altLoc - alternative location indicator. By default 1028 // (""), maximum-occupancy locations are taken. 1029 PChain chain; 1030 PPResidue res; 1031 PPAtom atom; 1032 int nResidues, nAtoms, i,j,k; 1033 1034 if (Calphas) { 1035 delete[] Calphas; 1036 Calphas = NULL; 1037 } 1038 nCalphas = 0; 1039 1040 if (residue) chain = residue->chain; 1041 else chain = NULL; 1042 1043 if (chain) { 1044 1045 chain->GetResidueTable ( res,nResidues ); 1046 1047 if (nResidues>0) { 1048 1049 Calphas = new PAtom[nResidues]; 1050 1051 if ((!altLoc[0]) || (altLoc[0]==' ')) { // main conformation 1052 1053 for (i=0;i<nResidues;i++) { 1054 nAtoms = res[i]->nAtoms; 1055 atom = res[i]->atom; 1056 for (j=0;j<nAtoms;j++) 1057 if (!strcmp(atom[j]->name," CA ")) { 1058 Calphas[nCalphas++] = atom[j]; 1059 break; 1060 } 1061 } 1062 1063 } else { // specific conformation 1064 1065 for (i=0;i<nResidues;i++) { 1066 nAtoms = res[i]->nAtoms; 1067 atom = res[i]->atom; 1068 k = 0; 1069 for (j=0;j<nAtoms;j++) 1070 if (!strcmp(atom[j]->name," CA ")) { 1071 k = 1; 1072 if (!atom[j]->altLoc[0]) { 1073 // take main conformation now as we do not know if 1074 // the specific one is in the file 1075 Calphas[nCalphas] = atom[j]; 1076 k = 2; 1077 } else if (atom[j]->altLoc[0]==altLoc[0]) { 1078 // get specific conformation and quit 1079 Calphas[nCalphas] = atom[j]; 1080 k = 2; 1081 break; 1082 } 1083 } else if (k) 1084 break; 1085 if (k==2) nCalphas++; 1086 } 1087 1088 } 1089 1090 } 1091 1092 } else if (residue) { // check just this atom's residue 1093 1094 Calphas = new PAtom[1]; // can't be more than 1 C-alpha! 1095 1096 nAtoms = residue->nAtoms; 1097 atom = residue->atom; 1098 1099 if ((!altLoc[0]) || (altLoc[0]==' ')) { // main conformation 1100 1101 for (j=0;j<nAtoms;j++) 1102 if (!strcmp(atom[j]->name," CA ")) { 1103 Calphas[nCalphas++] = atom[j]; 1104 break; 1105 } 1106 1107 } else { // specific conformation 1108 1109 k = 0; 1110 for (j=0;j<nAtoms;j++) 1111 if (!strcmp(atom[j]->name," CA ")) { 1112 k = 1; 1113 if (!atom[j]->altLoc[0]) { 1114 Calphas[nCalphas] = atom[j]; 1115 k = 2; 1116 } else if (atom[j]->altLoc[0]==altLoc[0]) { 1117 Calphas[nCalphas] = atom[j]; 1118 k = 2; 1119 break; 1120 } 1121 } else if (k) 1122 break; 1123 if (k==2) nCalphas++; 1124 1125 } 1126 1127 } 1128 1129 if (Calphas && (!nCalphas)) { 1130 delete[] Calphas; 1131 Calphas = NULL; 1132 } 1133 1134 } 1135 isMetal()1136 bool Atom::isMetal() { 1137 return mmdb::isMetal ( element ); 1138 } 1139 isSolvent()1140 bool Atom::isSolvent() { 1141 if (residue) return residue->isSolvent(); 1142 return false; 1143 } 1144 isInSelection(int selHnd)1145 bool Atom::isInSelection ( int selHnd ) { 1146 PRoot manager = (PRoot)GetCoordHierarchy(); 1147 PMask mask; 1148 if (manager) { 1149 mask = manager->GetSelMask ( selHnd ); 1150 if (mask) return CheckMask ( mask ); 1151 } 1152 return false; 1153 } 1154 isNTerminus()1155 bool Atom::isNTerminus() { 1156 if (residue) return residue->isNTerminus(); 1157 return false; 1158 } 1159 isCTerminus()1160 bool Atom::isCTerminus() { 1161 if (residue) return residue->isCTerminus(); 1162 return false; 1163 } 1164 CalAtomStatistics(RAtomStat AS)1165 void Atom::CalAtomStatistics ( RAtomStat AS ) { 1166 // AS must be initialized. The function only accumulates 1167 // the statistics. 1168 1169 if (!Ter) { 1170 1171 AS.nAtoms++; 1172 1173 if (AS.WhatIsSet & WhatIsSet & ASET_Coordinates) { 1174 GetStat ( x,AS.xmin,AS.xmax,AS.xm,AS.xm2 ); 1175 GetStat ( y,AS.ymin,AS.ymax,AS.ym,AS.ym2 ); 1176 GetStat ( z,AS.zmin,AS.zmax,AS.zm,AS.zm2 ); 1177 } else 1178 AS.WhatIsSet &= ~ASET_Coordinates; 1179 1180 if (AS.WhatIsSet & WhatIsSet & ASET_Occupancy) 1181 GetStat(occupancy,AS.occ_min,AS.occ_max,AS.occ_m,AS.occ_m2); 1182 else AS.WhatIsSet &= ~ASET_Occupancy; 1183 1184 if (AS.WhatIsSet & WhatIsSet & ASET_tempFactor) 1185 GetStat ( tempFactor,AS.tFmin,AS.tFmax,AS.tFm,AS.tFm2 ); 1186 else AS.WhatIsSet &= ~ASET_tempFactor; 1187 1188 if (AS.WhatIsSet & WhatIsSet & ASET_Anis_tFac) { 1189 GetStat ( u11,AS.u11_min,AS.u11_max,AS.u11_m,AS.u11_m2 ); 1190 GetStat ( u22,AS.u22_min,AS.u22_max,AS.u22_m,AS.u22_m2 ); 1191 GetStat ( u33,AS.u33_min,AS.u33_max,AS.u33_m,AS.u33_m2 ); 1192 GetStat ( u12,AS.u12_min,AS.u12_max,AS.u12_m,AS.u12_m2 ); 1193 GetStat ( u13,AS.u13_min,AS.u13_max,AS.u13_m,AS.u13_m2 ); 1194 GetStat ( u23,AS.u23_min,AS.u23_max,AS.u23_m,AS.u23_m2 ); 1195 } else 1196 AS.WhatIsSet &= ~ASET_Anis_tFac; 1197 1198 } 1199 1200 } 1201 1202 GetDist2(PAtom a)1203 realtype Atom::GetDist2 ( PAtom a ) { 1204 realtype dx,dy,dz; 1205 dx = a->x - x; 1206 dy = a->y - y; 1207 dz = a->z - z; 1208 return dx*dx + dy*dy + dz*dz; 1209 } 1210 GetDist2(PAtom a,mat44 & tm)1211 realtype Atom::GetDist2 ( PAtom a, mat44 & tm ) { 1212 realtype dx,dy,dz; 1213 dx = tm[0][0]*a->x + tm[0][1]*a->y + tm[0][2]*a->z + tm[0][3] - x; 1214 dy = tm[1][0]*a->x + tm[1][1]*a->y + tm[1][2]*a->z + tm[1][3] - y; 1215 dz = tm[2][0]*a->x + tm[2][1]*a->y + tm[2][2]*a->z + tm[2][3] - z; 1216 return dx*dx + dy*dy + dz*dz; 1217 } 1218 GetDist2(PAtom a,mat33 & r,vect3 & t)1219 realtype Atom::GetDist2 ( PAtom a, mat33 & r, vect3 & t ) { 1220 realtype dx,dy,dz; 1221 dx = r[0][0]*a->x + r[0][1]*a->y + r[0][2]*a->z + t[0] - x; 1222 dy = r[1][0]*a->x + r[1][1]*a->y + r[1][2]*a->z + t[1] - y; 1223 dz = r[2][0]*a->x + r[2][1]*a->y + r[2][2]*a->z + t[2] - z; 1224 return dx*dx + dy*dy + dz*dz; 1225 } 1226 GetDist2(realtype ax,realtype ay,realtype az)1227 realtype Atom::GetDist2 ( realtype ax, realtype ay, realtype az ) { 1228 realtype dx,dy,dz; 1229 dx = ax - x; 1230 dy = ay - y; 1231 dz = az - z; 1232 return dx*dx + dy*dy + dz*dz; 1233 } 1234 GetDist2(mat44 & tm,realtype ax,realtype ay,realtype az)1235 realtype Atom::GetDist2 ( mat44 & tm, // applies to 'this' 1236 realtype ax, realtype ay, realtype az ) { 1237 realtype dx,dy,dz; 1238 dx = tm[0][0]*x + tm[0][1]*y + tm[0][2]*z + tm[0][3] - ax; 1239 dy = tm[1][0]*x + tm[1][1]*y + tm[1][2]*z + tm[1][3] - ay; 1240 dz = tm[2][0]*x + tm[2][1]*y + tm[2][2]*z + tm[2][3] - az; 1241 return dx*dx + dy*dy + dz*dz; 1242 } 1243 GetDist2(vect3 & xyz)1244 realtype Atom::GetDist2 ( vect3 & xyz ) { 1245 realtype dx,dy,dz; 1246 dx = xyz[0] - x; 1247 dy = xyz[1] - y; 1248 dz = xyz[2] - z; 1249 return dx*dx + dy*dy + dz*dz; 1250 } 1251 GetCosine(PAtom a1,PAtom a2)1252 realtype Atom::GetCosine ( PAtom a1, PAtom a2 ) { 1253 // Calculates cosing of angle a1-this-a2, i.e. that between 1254 // bond [a1,this] and [this,a2]. 1255 realtype dx1,dy1,dz1, dx2,dy2,dz2,r; 1256 1257 dx1 = a1->x - x; 1258 dy1 = a1->y - y; 1259 dz1 = a1->z - z; 1260 r = dx1*dx1 + dy1*dy1 + dz1*dz1; 1261 1262 dx2 = a2->x - x; 1263 dy2 = a2->y - y; 1264 dz2 = a2->z - z; 1265 r *= dx2*dx2 + dy2*dy2 + dz2*dz2; 1266 1267 if (r>0.0) return (dx1*dx2 + dy1*dy2 + dz1*dz2)/sqrt(r); 1268 else return 0.0; 1269 1270 } 1271 1272 MakeTer()1273 void Atom::MakeTer() { 1274 WhatIsSet |= ASET_Coordinates; 1275 Het = false; 1276 Ter = true; 1277 } 1278 1279 SetAtomName(const AtomName atomName)1280 void Atom::SetAtomName ( const AtomName atomName ) { 1281 strcpy ( name,atomName ); 1282 } 1283 1284 SetElementName(const Element elName)1285 void Atom::SetElementName ( const Element elName ) { 1286 strcpy ( element,elName ); 1287 if (!element[0]) strcpy ( element," " ); 1288 else if ((!element[1]) || (element[1]==' ')) { 1289 element[2] = char(0); 1290 element[1] = element[0]; 1291 element[0] = ' '; 1292 } 1293 } 1294 SetCharge(cpstr chrg)1295 void Atom::SetCharge ( cpstr chrg ) { 1296 pstr p; 1297 charge = strtod ( chrg,&p ); 1298 if (p!=chrg) { 1299 WhatIsSet |= ASET_Charge; 1300 if ((charge>0.0) && (*p=='-')) 1301 charge = -charge; 1302 } 1303 } 1304 SetCharge(realtype chrg)1305 void Atom::SetCharge ( realtype chrg ) { 1306 if (chrg<MaxReal) { 1307 charge = chrg; 1308 WhatIsSet |= ASET_Charge; 1309 } 1310 } 1311 SetAtomIndex(int ix)1312 void Atom::SetAtomIndex ( int ix ) { 1313 // don't use in your applications! 1314 index = ix; 1315 } 1316 GetAtomID(pstr AtomID)1317 pstr Atom::GetAtomID ( pstr AtomID ) { 1318 char S[50]; 1319 AtomID[0] = char(0); 1320 if (residue) { 1321 if (residue->chain) { 1322 if (residue->chain->model) 1323 sprintf (AtomID,"/%i/",residue->chain->model->GetSerNum()); 1324 else strcpy ( AtomID,"/-/" ); 1325 strcat ( AtomID,residue->chain->chainID ); 1326 } else 1327 strcpy ( AtomID,"/-/-" ); 1328 ParamStr ( AtomID,pstr("/"),residue->seqNum ); 1329 if (residue->name[0]) { 1330 strcat ( AtomID,"(" ); 1331 strcat ( AtomID,residue->name ); 1332 strcat ( AtomID,")" ); 1333 } 1334 if (residue->insCode[0]) { 1335 strcat ( AtomID,"." ); 1336 strcat ( AtomID,residue->insCode ); 1337 } 1338 strcat ( AtomID,"/" ); 1339 } else 1340 strcpy ( AtomID,"/-/-/-/" ); 1341 strcpy_css ( S,name ); 1342 if (!S[0]) strcpy ( S,"-" ); 1343 strcat ( AtomID,S ); 1344 strcpy_css ( S,element ); 1345 if (S[0]) { 1346 strcat ( AtomID,"[" ); 1347 strcat ( AtomID,S ); 1348 strcat ( AtomID,"]" ); 1349 } 1350 if (altLoc[0]) { 1351 strcat ( AtomID,":" ); 1352 strcat ( AtomID,altLoc ); 1353 } 1354 return AtomID; 1355 } 1356 GetAtomIDfmt(pstr AtomID)1357 pstr Atom::GetAtomIDfmt ( pstr AtomID ) { 1358 int n; 1359 char S[50]; 1360 AtomID[0] = char(0); 1361 if (residue) { 1362 if (residue->chain) { 1363 if (residue->chain->model) { 1364 n = residue->chain->model->GetNumberOfModels(); 1365 if (n<10) strcpy ( S,"/%1i/" ); 1366 else if (n<100) strcpy ( S,"/%2i/" ); 1367 else if (n<1000) strcpy ( S,"/%3i/" ); 1368 else strcpy ( S,"/%i/" ); 1369 sprintf ( AtomID,S,residue->chain->model->GetSerNum() ); 1370 } else 1371 strcpy ( AtomID,"/-/" ); 1372 strcat ( AtomID,residue->chain->chainID ); 1373 } else 1374 strcpy ( AtomID,"/-/-" ); 1375 if ((-999<=residue->seqNum) && (residue->seqNum<=9999)) 1376 sprintf ( S,"/%4i",residue->seqNum ); 1377 else sprintf ( S,"/%i" ,residue->seqNum ); 1378 strcat ( AtomID,S ); 1379 sprintf ( S,"(%3s).%1s/",residue->name,residue->insCode ); 1380 strcat ( AtomID,S ); 1381 } else 1382 strcpy ( AtomID,"/-/-/----(---).-/" ); 1383 sprintf ( S,"%4s[%2s]:%1s",name,element,altLoc ); 1384 strcat ( AtomID,S ); 1385 return AtomID; 1386 } 1387 1388 1389 ConvertPDBHETATM(int ix,cpstr S)1390 ERROR_CODE Atom::ConvertPDBHETATM ( int ix, cpstr S ) { 1391 // Gets data from the PDB ASCII HETATM record. 1392 // This function DOES NOT check the "HETATM" keyword and 1393 // does not decode the chain and residue parameters! These 1394 // must be treated by the calling process, see 1395 // Chain::ConvertPDBASCII(). 1396 ERROR_CODE RC; 1397 RC = ConvertPDBATOM ( ix,S ); 1398 Het = true; 1399 return RC; 1400 } 1401 GetData(cpstr S)1402 void Atom::GetData ( cpstr S ) { 1403 pstr p; 1404 1405 if (((S[6]>='0') && (S[6]<='9')) || (S[6]==' ')) { 1406 // Here we forgive cards with unreadable serial numbers 1407 // as we always have index (ix) for the card. For the sake 1408 // of strict PDB syntax we would have to return 1409 // Error_UnrecognizedInteger here. 1410 if (!(GetInteger(serNum,&(S[6]),5))) serNum = -1; 1411 } else 1412 hy36decode ( 5,&(S[6]),5,&serNum ); 1413 1414 // if (!(GetInteger(serNum,&(S[6]),5))) serNum = index; 1415 1416 altLoc[0] = S[16]; 1417 if (altLoc[0]==' ') altLoc[0] = char(0); 1418 else altLoc[1] = char(0); 1419 GetString ( name ,&(S[12]),4 ); 1420 strcpy_ncss ( segID ,&(S[72]),4 ); 1421 GetString ( element,&(S[76]),2 ); 1422 charge = strtod ( &(S[78]),&p ); 1423 if ((charge!=0.0) && (p!=&(S[78]))) { 1424 WhatIsSet |= ASET_Charge; 1425 if ((charge>0.0) && (*p=='-')) 1426 charge = -charge; 1427 } 1428 1429 RestoreElementName(); 1430 strcpy ( label_atom_id,name ); 1431 1432 } 1433 CheckData(cpstr S)1434 ERROR_CODE Atom::CheckData ( cpstr S ) { 1435 int sN; 1436 AltLoc aloc; 1437 SegID sID; 1438 Element elmnt; 1439 pstr p; 1440 realtype achrg; 1441 1442 aloc[0] = S[16]; 1443 if (aloc[0]==' ') aloc[0] = char(0); 1444 else aloc[1] = char(0); 1445 1446 strcpy_ncss ( sID ,&(S[72]),4 ); 1447 GetString ( elmnt,&(S[76]),2 ); 1448 1449 if (ignoreCharge) 1450 achrg = charge; 1451 else { 1452 achrg = strtod ( &(S[78]),&p ); 1453 if ((achrg!=0.0) && (p!=&(S[78]))) { 1454 if ((achrg>0.0) && (*p=='-')) 1455 achrg = -achrg; 1456 } 1457 } 1458 1459 // if (!(GetInteger(sN,&(S[6]),5))) 1460 // sN = index; 1461 1462 if (hy36decode(5,&(S[6]),5,&sN)) 1463 sN = index; 1464 1465 if (ignoreSegID) { 1466 if (segID[0]) strcpy ( sID,segID ); 1467 else strcpy ( segID,sID ); 1468 } 1469 1470 if (ignoreElement) { 1471 if (element[0]) strcpy ( elmnt,element ); 1472 else strcpy ( element,elmnt ); 1473 } 1474 1475 if (ignoreUnmatch) return Error_NoError; 1476 1477 // Here we forgive cards with unreadable serial numbers 1478 // as we always have index (ix) for the card. For the sake 1479 // of strict PDB syntax we would have to return 1480 // Error_UnrecognizedInteger . 1481 if ((sN!=serNum) || 1482 (strcmp (altLoc ,aloc )) || 1483 (strncmp(name ,&(S[12]),4)) || 1484 (strcmp (segID ,sID )) || 1485 (strcmp (element,elmnt )) || 1486 (charge!=achrg)) { 1487 /* 1488 char name1[100]; 1489 strncpy ( name1,&(S[12]),4 ); name1[4] = char(0); 1490 printf ( "\n serNum %5i %5i\n" 1491 " residue '%s' '%s'\n" 1492 " altLoc '%s' '%s'\n" 1493 " name '%s' '%s'\n" 1494 " segId '%s' '%s'\n" 1495 " element '%s' '%s'\n" 1496 " charge '%s' '%s'\n", 1497 sN,serNum, res->name,residue->name, 1498 altLoc ,aloc, name,name1, 1499 segID ,sID, 1500 element,elmnt, 1501 charge ,achrg ); 1502 if (res!=residue) printf (" it's a residue\n" ); 1503 */ 1504 return Error_ATOM_Unmatch; 1505 } 1506 1507 return Error_NoError; 1508 1509 } 1510 1511 GetCIF(int ix,mmcif::PLoop Loop,mmcif::PLoop LoopAnis)1512 ERROR_CODE Atom::GetCIF ( int ix, mmcif::PLoop Loop, 1513 mmcif::PLoop LoopAnis ) { 1514 char PDBGroup[30]; 1515 int k; 1516 ERROR_CODE RC; 1517 1518 index = ix; 1519 1520 if (WhatIsSet & ASET_Coordinates) 1521 return Error_ATOM_AlreadySet; 1522 1523 /* 1524 1525 loop_ 1526 *0 _atom_site.group_PDB 1527 *1 _atom_site.id 1528 *2 _atom_site.type_symbol <- chem elem 1529 -3 _atom_site.label_atom_id <- atom name 1530 *4 _atom_site.label_alt_id <- alt code 1531 =5 _atom_site.label_comp_id <- res name ??? 1532 =6 _atom_site.label_asym_id <- chain id ??? 1533 =7 _atom_site.label_entity_id < ??? 1534 =8 _atom_site.label_seq_id <- poly seq id 1535 +9 _atom_site.pdbx_PDB_ins_code <- ins code 1536 -10 _atom_site.segment_id <- segment id 1537 *11 _atom_site.Cartn_x 1538 *12 _atom_site.Cartn_y 1539 *13 _atom_site.Cartn_z 1540 *14 _atom_site.occupancy 1541 *15 _atom_site.B_iso_or_equiv 1542 *16 _atom_site.Cartn_x_esd 1543 *17 _atom_site.Cartn_y_esd 1544 *18 _atom_site.Cartn_z_esd 1545 *19 _atom_site.occupancy_esd 1546 *20 _atom_site.B_iso_or_equiv_esd 1547 *21 _atom_site.pdbx_formal_charge 1548 +22 _atom_site.auth_seq_id <- seq id we want 1549 +23 _atom_site.auth_comp_id <- res name we want 1550 +24 _atom_site.auth_asym_id <- ch id we want ? 1551 *25 _atom_site.auth_atom_id <- atom name we want ? 1552 +26 _atom_site.pdbx_PDB_model_num <- model no 1553 1554 '+' read in Root::CheckAtomPlace() 1555 '=' new in residue, read in Root::CheckAtomPlace() 1556 '-' new in atom 1557 1558 */ 1559 1560 1561 // (0) 1562 k = ix-1; 1563 CIFGetString ( PDBGroup,Loop,CIFTAG_GROUP_PDB,k, 1564 sizeof(PDBGroup),pstr("") ); 1565 1566 Ter = !strcmp(PDBGroup,pstr("TER") ); 1567 Het = !strcmp(PDBGroup,pstr("HETATM")); 1568 1569 // (1) 1570 RC = CIFGetInteger1 ( serNum,Loop,CIFTAG_ID,k ); 1571 if (RC) { 1572 if (Ter) serNum = -1; 1573 else if (RC==Error_NoData) serNum = index; 1574 else 1575 return RC; 1576 } 1577 1578 if (Ter) { 1579 Loop->DeleteRow ( k ); 1580 WhatIsSet |= ASET_Coordinates; 1581 return Error_NoError; 1582 } 1583 1584 // (25) 1585 CIFGetString ( name,Loop,CIFTAG_AUTH_ATOM_ID,k, 1586 sizeof(name) ,pstr("") ); 1587 // (3) 1588 CIFGetString ( label_atom_id,Loop,CIFTAG_LABEL_ATOM_ID,k, 1589 sizeof(label_atom_id),pstr("") ); 1590 if (!name[0]) 1591 strcpy ( name,label_atom_id ); 1592 // (4) 1593 CIFGetString ( altLoc,Loop,CIFTAG_LABEL_ALT_ID ,k, 1594 sizeof(altLoc),pstr("") ); 1595 1596 // (11,12,13) 1597 RC = CIFGetReal1 ( x,Loop,CIFTAG_CARTN_X,k ); 1598 if (!RC) RC = CIFGetReal1 ( y,Loop,CIFTAG_CARTN_Y,k ); 1599 if (!RC) RC = CIFGetReal1 ( z,Loop,CIFTAG_CARTN_Z,k ); 1600 if (RC) return Error_ATOM_Unrecognized; 1601 WhatIsSet |= ASET_Coordinates; 1602 1603 // (14) 1604 if (!CIFGetReal1(occupancy,Loop,CIFTAG_OCCUPANCY,k)) 1605 WhatIsSet |= ASET_Occupancy; 1606 // (15) 1607 if (!CIFGetReal1(tempFactor,Loop,CIFTAG_B_ISO_OR_EQUIV,k)) 1608 WhatIsSet |= ASET_tempFactor; 1609 1610 // (10) 1611 CIFGetString ( segID,Loop,CIFTAG_SEGMENT_ID,k, 1612 sizeof(segID) ,pstr("") ); 1613 // (21) 1614 if (!CIFGetReal1(charge,Loop,CIFTAG_PDBX_FORMAL_CHARGE,k)) 1615 WhatIsSet |= ASET_Charge; 1616 // (2) 1617 RC = CIFGetString ( element,Loop,CIFTAG_TYPE_SYMBOL,k, 1618 sizeof(element),pstr("") ); 1619 if (RC) 1620 CIFGetString ( element,Loop,CIFTAG_ATOM_TYPE_SYMBOL,k, 1621 sizeof(element),pstr("") ); 1622 1623 RestoreElementName(); 1624 MakePDBAtomName(); 1625 1626 // (16,17,18) 1627 RC = CIFGetReal1 ( sigX,Loop,CIFTAG_CARTN_X_ESD,k ); 1628 if (!RC) RC = CIFGetReal1 ( sigY,Loop,CIFTAG_CARTN_Y_ESD,k ); 1629 if (!RC) RC = CIFGetReal1 ( sigZ,Loop,CIFTAG_CARTN_Z_ESD,k ); 1630 if (RC==Error_UnrecognizedReal) 1631 return RC; 1632 if (!RC) WhatIsSet |= ASET_CoordSigma; 1633 1634 // (19) 1635 if (!CIFGetReal1(sigOcc,Loop,CIFTAG_OCCUPANCY_ESD,k)) 1636 WhatIsSet |= ASET_OccSigma; 1637 // (20) 1638 if (!CIFGetReal1(sigTemp,Loop,CIFTAG_B_ISO_OR_EQUIV_ESD,k)) 1639 WhatIsSet |= ASET_tFacSigma; 1640 1641 Loop->DeleteRow ( k ); 1642 1643 if (LoopAnis) { 1644 1645 RC = CIFGetReal1 ( u11,LoopAnis,CIFTAG_U11,k ); 1646 if (!RC) RC = CIFGetReal1 ( u22,LoopAnis,CIFTAG_U22,k ); 1647 if (!RC) RC = CIFGetReal1 ( u33,LoopAnis,CIFTAG_U33,k ); 1648 if (!RC) RC = CIFGetReal1 ( u13,LoopAnis,CIFTAG_U13,k ); 1649 if (!RC) RC = CIFGetReal1 ( u12,LoopAnis,CIFTAG_U12,k ); 1650 if (!RC) RC = CIFGetReal1 ( u23,LoopAnis,CIFTAG_U23,k ); 1651 if (RC==Error_UnrecognizedReal) 1652 return RC; 1653 if (!RC) WhatIsSet |= ASET_Anis_tFac; 1654 1655 RC = CIFGetReal1 ( su11,LoopAnis,CIFTAG_U11_ESD,k ); 1656 if (!RC) RC = CIFGetReal1 ( su22,LoopAnis,CIFTAG_U22_ESD,k ); 1657 if (!RC) RC = CIFGetReal1 ( su33,LoopAnis,CIFTAG_U33_ESD,k ); 1658 if (!RC) RC = CIFGetReal1 ( su13,LoopAnis,CIFTAG_U13_ESD,k ); 1659 if (!RC) RC = CIFGetReal1 ( su12,LoopAnis,CIFTAG_U12_ESD,k ); 1660 if (!RC) RC = CIFGetReal1 ( su23,LoopAnis,CIFTAG_U23_ESD,k ); 1661 if (RC==Error_UnrecognizedReal) 1662 return RC; 1663 if (!RC) WhatIsSet |= ASET_Anis_tFSigma; 1664 1665 LoopAnis->DeleteRow ( k ); 1666 1667 } 1668 1669 return Error_NoError; 1670 1671 } 1672 RestoreElementName()1673 bool Atom::RestoreElementName() { 1674 // This function works only if element name is not given. 1675 1676 if (Ter) { 1677 name[0] = char(0); 1678 element[0] = char(0); 1679 return false; 1680 } 1681 if ((!element[0]) || 1682 ((element[0]==' ') && ((!element[1]) || (element[1]==' ')))) { 1683 if (strlen(name)==4) { 1684 if ((name[0]>='A') && (name[0]<='Z')) { 1685 element[0] = name[0]; 1686 element[1] = name[1]; 1687 } else { 1688 element[0] = ' '; 1689 element[1] = name[1]; 1690 } 1691 } else { // nasty hack for mmcif files without element column 1692 element[0] = ' '; 1693 element[1] = name[0]; 1694 } 1695 element[2] = char(0); 1696 return false; 1697 } else if (!element[1]) { 1698 // not aligned element name, possibly coming from mmCIF 1699 element[1] = element[0]; 1700 element[0] = ' '; 1701 element[2] = char(0); 1702 return false; 1703 } 1704 1705 return true; 1706 1707 } 1708 1709 MakePDBAtomName()1710 bool Atom::MakePDBAtomName() { 1711 int i,k; 1712 1713 if (Ter) { 1714 name [0] = char(0); 1715 element[0] = char(0); 1716 return false; 1717 } 1718 1719 UpperCase ( name ); 1720 UpperCase ( element ); 1721 1722 if (strlen(name)>=4) 1723 return true; // nothing to do 1724 1725 if ((element[0]==' ') && (element[1]==' ')) { 1726 // element name not given, make one from the atom name 1727 if ((name[0]>='A') && (name[0]<='Z')) { 1728 if (!name[1]) { 1729 name[4] = char(0); 1730 name[3] = ' '; 1731 name[2] = ' '; 1732 name[1] = name[0]; 1733 name[0] = ' '; 1734 } 1735 /* the commented part looks like a wrong inheritance 1736 from FORTRAN RWBrook. Commented on 04.03.2004, 1737 to be removed. 1738 else if ((name[0]=='C') && (name[1]=='A')) { 1739 name[4] = char(0); 1740 name[3] = name[2]; 1741 name[2] = name[1]; 1742 name[1] = name[0]; 1743 name[0] = ' '; 1744 } 1745 */ 1746 element[0] = name[0]; 1747 element[1] = name[1]; 1748 } else { 1749 element[0] = ' '; 1750 element[1] = name[1]; 1751 } 1752 element[2] = char(0); 1753 return false; 1754 } else if ((name[0]>='A') && (name[0]<='Z')) { 1755 if (!element[1]) { 1756 element[2] = char(0); 1757 element[1] = element[0]; 1758 element[0] = ' '; 1759 k = strlen(name); 1760 if (k<4) { 1761 for (i=3;i>0;i--) 1762 name[i] = name[i-1]; 1763 name[0] = ' '; 1764 k++; 1765 while (k<4) 1766 name[k++] = ' '; 1767 name[k] = char(0); 1768 } 1769 } else if ((element[0]==' ') && (element[1]!=name[1])) { 1770 for (i=3;i>0;i--) 1771 name[i] = name[i-1]; 1772 name[0] = ' '; 1773 name[4] = char(0); 1774 k = strlen(name); 1775 while (k<4) 1776 name[k++] = ' '; 1777 name[k] = char(0); 1778 } else { 1779 k = strlen(name); 1780 while (k<4) 1781 name[k++] = ' '; 1782 name[k] = char(0); 1783 } 1784 } 1785 return true; 1786 } 1787 SetCoordinates(realtype xx,realtype yy,realtype zz,realtype occ,realtype tFac)1788 void Atom::SetCoordinates ( realtype xx, realtype yy, realtype zz, 1789 realtype occ, realtype tFac ) { 1790 x = xx; 1791 y = yy; 1792 z = zz; 1793 occupancy = occ; 1794 tempFactor = tFac; 1795 WhatIsSet |= ASET_Coordinates | ASET_Occupancy | ASET_tempFactor; 1796 } 1797 Transform(const mat33 & tm,vect3 & v)1798 void Atom::Transform ( const mat33 & tm, vect3 & v ) { 1799 realtype x1,y1,z1; 1800 x1 = tm[0][0]*x + tm[0][1]*y + tm[0][2]*z + v[0]; 1801 y1 = tm[1][0]*x + tm[1][1]*y + tm[1][2]*z + v[1]; 1802 z1 = tm[2][0]*x + tm[2][1]*y + tm[2][2]*z + v[2]; 1803 x = x1; 1804 y = y1; 1805 z = z1; 1806 } 1807 Transform(const mat44 & tm)1808 void Atom::Transform ( const mat44 & tm ) { 1809 realtype x1,y1,z1; 1810 x1 = tm[0][0]*x + tm[0][1]*y + tm[0][2]*z + tm[0][3]; 1811 y1 = tm[1][0]*x + tm[1][1]*y + tm[1][2]*z + tm[1][3]; 1812 z1 = tm[2][0]*x + tm[2][1]*y + tm[2][2]*z + tm[2][3]; 1813 x = x1; 1814 y = y1; 1815 z = z1; 1816 } 1817 TransformCopy(const mat44 & tm,realtype & xx,realtype & yy,realtype & zz)1818 void Atom::TransformCopy ( const mat44 & tm, 1819 realtype & xx, 1820 realtype & yy, 1821 realtype & zz ) { 1822 xx = tm[0][0]*x + tm[0][1]*y + tm[0][2]*z + tm[0][3]; 1823 yy = tm[1][0]*x + tm[1][1]*y + tm[1][2]*z + tm[1][3]; 1824 zz = tm[2][0]*x + tm[2][1]*y + tm[2][2]*z + tm[2][3]; 1825 } 1826 TransformCopy(const mat44 & tm,vect3 & xyz)1827 void Atom::TransformCopy ( const mat44 & tm, vect3 & xyz ) { 1828 xyz[0] = tm[0][0]*x + tm[0][1]*y + tm[0][2]*z + tm[0][3]; 1829 xyz[1] = tm[1][0]*x + tm[1][1]*y + tm[1][2]*z + tm[1][3]; 1830 xyz[2] = tm[2][0]*x + tm[2][1]*y + tm[2][2]*z + tm[2][3]; 1831 } 1832 TransformSet(const mat44 & tm,realtype xx,realtype yy,realtype zz)1833 void Atom::TransformSet ( const mat44 & tm, 1834 realtype xx, 1835 realtype yy, 1836 realtype zz ) { 1837 x = tm[0][0]*xx + tm[0][1]*yy + tm[0][2]*zz + tm[0][3]; 1838 y = tm[1][0]*xx + tm[1][1]*yy + tm[1][2]*zz + tm[1][3]; 1839 z = tm[2][0]*xx + tm[2][1]*yy + tm[2][2]*zz + tm[2][3]; 1840 } 1841 1842 1843 // ------- user-defined data handlers 1844 PutUDData(int UDDhandle,int iudd)1845 int Atom::PutUDData ( int UDDhandle, int iudd ) { 1846 if (UDDhandle & UDRF_ATOM) 1847 return UDData::putUDData ( UDDhandle,iudd ); 1848 else return UDDATA_WrongUDRType; 1849 } 1850 PutUDData(int UDDhandle,realtype rudd)1851 int Atom::PutUDData ( int UDDhandle, realtype rudd ) { 1852 if (UDDhandle & UDRF_ATOM) 1853 return UDData::putUDData ( UDDhandle,rudd ); 1854 else return UDDATA_WrongUDRType; 1855 } 1856 PutUDData(int UDDhandle,cpstr sudd)1857 int Atom::PutUDData ( int UDDhandle, cpstr sudd ) { 1858 if (UDDhandle & UDRF_ATOM) 1859 return UDData::putUDData ( UDDhandle,sudd ); 1860 else return UDDATA_WrongUDRType; 1861 } 1862 GetUDData(int UDDhandle,int & iudd)1863 int Atom::GetUDData ( int UDDhandle, int & iudd ) { 1864 if (UDDhandle & UDRF_ATOM) 1865 return UDData::getUDData ( UDDhandle,iudd ); 1866 else return UDDATA_WrongUDRType; 1867 } 1868 GetUDData(int UDDhandle,realtype & rudd)1869 int Atom::GetUDData ( int UDDhandle, realtype & rudd ) { 1870 if (UDDhandle & UDRF_ATOM) 1871 return UDData::getUDData ( UDDhandle,rudd ); 1872 else return UDDATA_WrongUDRType; 1873 } 1874 GetUDData(int UDDhandle,pstr sudd,int maxLen)1875 int Atom::GetUDData ( int UDDhandle, pstr sudd, int maxLen ) { 1876 if (UDDhandle & UDRF_ATOM) 1877 return UDData::getUDData ( UDDhandle,sudd,maxLen ); 1878 else return UDDATA_WrongUDRType; 1879 } 1880 GetUDData(int UDDhandle,pstr & sudd)1881 int Atom::GetUDData ( int UDDhandle, pstr & sudd ) { 1882 if (UDDhandle & UDRF_ATOM) 1883 return UDData::getUDData ( UDDhandle,sudd ); 1884 else return UDDATA_WrongUDRType; 1885 } 1886 1887 Copy(PAtom atom)1888 void Atom::Copy ( PAtom atom ) { 1889 // this does not make any references in residues and does 1890 // not change indices!! it does change serial numbers, though. 1891 1892 serNum = atom->serNum; 1893 x = atom->x; 1894 y = atom->y; 1895 z = atom->z; 1896 occupancy = atom->occupancy; 1897 tempFactor = atom->tempFactor; 1898 sigX = atom->sigX; 1899 sigY = atom->sigY; 1900 sigZ = atom->sigZ; 1901 sigOcc = atom->sigOcc; 1902 sigTemp = atom->sigTemp; 1903 u11 = atom->u11; 1904 u22 = atom->u22; 1905 u33 = atom->u33; 1906 u12 = atom->u12; 1907 u13 = atom->u13; 1908 u23 = atom->u23; 1909 su11 = atom->su11; 1910 su22 = atom->su22; 1911 su33 = atom->su33; 1912 su12 = atom->su12; 1913 su13 = atom->su13; 1914 su23 = atom->su23; 1915 Het = atom->Het; 1916 Ter = atom->Ter; 1917 WhatIsSet = atom->WhatIsSet; 1918 1919 strcpy ( name ,atom->name ); 1920 strcpy ( label_atom_id,atom->label_atom_id ); 1921 strcpy ( altLoc ,atom->altLoc ); 1922 strcpy ( segID ,atom->segID ); 1923 strcpy ( element ,atom->element ); 1924 strcpy ( energyType ,atom->energyType ); 1925 charge = atom->charge; 1926 1927 } 1928 CheckID(const AtomName aname,const Element elname,const AltLoc aloc)1929 int Atom::CheckID ( const AtomName aname, const Element elname, 1930 const AltLoc aloc ) { 1931 pstr p1,p2; 1932 if (aname) { 1933 if (aname[0]!='*') { 1934 p1 = name; 1935 while (*p1==' ') p1++; 1936 p2 = pstr(aname); 1937 while (*p2==' ') p2++; 1938 while ((*p2) && (*p1) && (*p1!=' ') && (*p2!=' ')) { 1939 if (*p1!=*p2) return 0; 1940 p1++; 1941 p2++; 1942 } 1943 if (*p1!=*p2) { 1944 if (((*p1) && (*p1!=' ')) || 1945 ((*p2) && (*p2!=' '))) return 0; 1946 } 1947 } 1948 } 1949 if (elname) { 1950 if (elname[0]!='*') { 1951 p1 = element; 1952 while (*p1==' ') p1++; 1953 p2 = pstr(elname); 1954 while (*p2==' ') p2++; 1955 while ((*p2) && (*p1) && (*p1!=' ') && (*p2!=' ')) { 1956 if (*p1!=*p2) return 0; 1957 p1++; 1958 p2++; 1959 } 1960 if (*p1!=*p2) return 0; 1961 } 1962 } 1963 if (aloc) { 1964 if ((aloc[0]!='*') && (strcmp(aloc,altLoc))) return 0; 1965 } 1966 return 1; 1967 } 1968 CheckIDS(cpstr ID)1969 int Atom::CheckIDS ( cpstr ID ) { 1970 AtomName aname; 1971 Element elname; 1972 AltLoc aloc; 1973 pstr p; 1974 p = LastOccurence ( ID,'/' ); 1975 if (p) p++; 1976 else p = pstr(ID); 1977 ParseAtomID ( p,aname,elname,aloc ); 1978 return CheckID ( aname,elname,aloc ); 1979 } 1980 SetCompactBinary()1981 void Atom::SetCompactBinary() { 1982 WhatIsSet |= ASET_CompactBinary; 1983 } 1984 write(io::RFile f)1985 void Atom::write ( io::RFile f ) { 1986 int i,k; 1987 byte Version=2; 1988 byte nb; 1989 1990 f.WriteWord ( &WhatIsSet ); 1991 if (WhatIsSet & ASET_CompactBinary) { 1992 if (Ter) WhatIsSet |= ASET_ShortTer; 1993 if (Het) WhatIsSet |= ASET_ShortHet; 1994 f.WriteInt ( &index ); 1995 f.WriteTerLine ( name ,false ); 1996 f.WriteTerLine ( altLoc ,false ); 1997 f.WriteTerLine ( element,false ); 1998 if (WhatIsSet & ASET_Coordinates) { 1999 /* 2000 char S[100]; 2001 sprintf ( S,"%8.3f%8.3f%8.3f",x,y,z ); 2002 f.WriteFile ( S,24 ); 2003 */ 2004 /* 2005 f.WriteReal ( &x ); 2006 f.WriteReal ( &y ); 2007 f.WriteReal ( &z ); 2008 */ 2009 2010 // Make integer conversion in order to eliminate round-offs 2011 // from PDB format %8.3f for coordinates and save space 2012 // comparing to real*8 numbers. Writing floats causes 2013 // precision loss 2014 k = mround(x*10000.0); f.WriteInt ( &k ); 2015 k = mround(y*10000.0); f.WriteInt ( &k ); 2016 k = mround(z*10000.0); f.WriteInt ( &k ); 2017 2018 /* 2019 f.WriteFloat ( &x ); 2020 f.WriteFloat ( &y ); 2021 f.WriteFloat ( &z ); 2022 */ 2023 } 2024 return; 2025 } 2026 2027 f.WriteByte ( &Version ); 2028 2029 UDData::write ( f ); 2030 2031 f.WriteInt ( &serNum ); 2032 f.WriteInt ( &index ); 2033 f.WriteTerLine ( name ,false ); 2034 f.WriteTerLine ( label_atom_id,false ); 2035 f.WriteTerLine ( altLoc ,false ); 2036 f.WriteTerLine ( segID ,false ); 2037 f.WriteTerLine ( element ,false ); 2038 f.WriteTerLine ( energyType ,false ); 2039 f.WriteFloat ( &charge ); 2040 f.WriteBool ( &Het ); 2041 f.WriteBool ( &Ter ); 2042 2043 if (WhatIsSet & ASET_Coordinates) { 2044 f.WriteFloat ( &x ); 2045 f.WriteFloat ( &y ); 2046 f.WriteFloat ( &z ); 2047 if (WhatIsSet & ASET_Occupancy) 2048 f.WriteFloat ( &occupancy ); 2049 if (WhatIsSet & ASET_tempFactor) 2050 f.WriteFloat ( &tempFactor ); 2051 } 2052 2053 if (WhatIsSet & ASET_CoordSigma) { 2054 f.WriteFloat ( &sigX ); 2055 f.WriteFloat ( &sigY ); 2056 f.WriteFloat ( &sigZ ); 2057 if ((WhatIsSet & ASET_Occupancy) && 2058 (WhatIsSet & ASET_OccSigma)) 2059 f.WriteFloat ( &sigOcc ); 2060 if ((WhatIsSet & ASET_tempFactor) && 2061 (WhatIsSet & ASET_tFacSigma)) 2062 f.WriteFloat ( &sigTemp ); 2063 } 2064 2065 if (WhatIsSet & ASET_Anis_tFac) { 2066 f.WriteFloat ( &u11 ); 2067 f.WriteFloat ( &u22 ); 2068 f.WriteFloat ( &u33 ); 2069 f.WriteFloat ( &u12 ); 2070 f.WriteFloat ( &u13 ); 2071 f.WriteFloat ( &u23 ); 2072 if (WhatIsSet & ASET_Anis_tFSigma) { 2073 f.WriteFloat ( &su11 ); 2074 f.WriteFloat ( &su22 ); 2075 f.WriteFloat ( &su33 ); 2076 f.WriteFloat ( &su12 ); 2077 f.WriteFloat ( &su13 ); 2078 f.WriteFloat ( &su23 ); 2079 } 2080 } 2081 2082 nb = byte(nBonds & 0x000000FF); 2083 f.WriteByte ( &nb ); 2084 for (i=0;i<nb;i++) 2085 if (Bond[i].atom) { 2086 f.WriteInt ( &(Bond[i].atom->index) ); 2087 f.WriteByte ( &(Bond[i].order) ); 2088 } else { 2089 k = -1; 2090 f.WriteInt ( &k ); 2091 } 2092 2093 } 2094 read(io::RFile f)2095 void Atom::read ( io::RFile f ) { 2096 int i,k; 2097 byte nb,Version; 2098 2099 FreeMemory(); 2100 2101 f.ReadWord ( &WhatIsSet ); 2102 if (WhatIsSet & ASET_CompactBinary) { 2103 f.ReadInt ( &index ); 2104 f.ReadTerLine ( name ,false ); 2105 f.ReadTerLine ( altLoc ,false ); 2106 f.ReadTerLine ( element,false ); 2107 if (WhatIsSet & ASET_Coordinates) { 2108 /* 2109 char S[100]; 2110 f.ReadFile ( S,24 ); 2111 GetReal(x,&(S[0]),8); 2112 GetReal(y,&(S[8]),8); 2113 GetReal(z,&(S[16]),8); 2114 */ 2115 /* 2116 f.ReadReal ( &x ); 2117 f.ReadReal( &y ); 2118 f.ReadReal ( &z ); 2119 */ 2120 2121 // Make integer conversion in order to eliminate round-offs 2122 // from PDB format %8.3f for coordinates and save space 2123 // comparing to real*8 numbers. Writing floats causes 2124 // precision loss 2125 f.ReadInt ( &k ); x = realtype(k)/10000.0; 2126 f.ReadInt ( &k ); y = realtype(k)/10000.0; 2127 f.ReadInt ( &k ); z = realtype(k)/10000.0; 2128 2129 /* 2130 f.ReadFloat ( &x ); 2131 f.ReadFloat ( &y ); 2132 f.ReadFloat ( &z ); 2133 */ 2134 } 2135 serNum = index; 2136 Ter = WhatIsSet & ASET_ShortTer; 2137 Het = WhatIsSet & ASET_ShortHet; 2138 name [4] = char(0); 2139 altLoc [1] = char(0); 2140 element[2] = char(0); 2141 segID [0] = char(0); 2142 charge = 0.0; 2143 WhatIsSet &= ASET_All; 2144 return; 2145 } 2146 2147 f.ReadByte ( &Version ); 2148 2149 UDData::read ( f ); 2150 2151 f.ReadInt ( &serNum ); 2152 f.ReadInt ( &index ); 2153 f.ReadTerLine ( name ,false ); 2154 if (Version>1) 2155 f.ReadTerLine ( label_atom_id,false ); 2156 f.ReadTerLine ( altLoc ,false ); 2157 f.ReadTerLine ( segID ,false ); 2158 f.ReadTerLine ( element ,false ); 2159 f.ReadTerLine ( energyType,false ); 2160 f.ReadFloat ( &charge ); 2161 f.ReadBool ( &Het ); 2162 f.ReadBool ( &Ter ); 2163 2164 if (WhatIsSet & ASET_Coordinates) { 2165 f.ReadFloat ( &x ); 2166 f.ReadFloat ( &y ); 2167 f.ReadFloat ( &z ); 2168 if (WhatIsSet & ASET_Occupancy) f.ReadFloat ( &occupancy ); 2169 else occupancy = 0.0; 2170 if (WhatIsSet & ASET_tempFactor) f.ReadFloat ( &tempFactor ); 2171 else tempFactor = 0.0; 2172 } else { 2173 x = 0.0; 2174 y = 0.0; 2175 z = 0.0; 2176 occupancy = 0.0; 2177 tempFactor = 0.0; 2178 } 2179 2180 if (WhatIsSet & ASET_CoordSigma) { 2181 f.ReadFloat ( &sigX ); 2182 f.ReadFloat ( &sigY ); 2183 f.ReadFloat ( &sigZ ); 2184 if ((WhatIsSet & ASET_Occupancy) && 2185 (WhatIsSet & ASET_OccSigma)) 2186 f.ReadFloat ( &sigOcc ); 2187 else sigOcc = 0.0; 2188 if ((WhatIsSet & ASET_tempFactor) && 2189 (WhatIsSet & ASET_tFacSigma)) 2190 f.ReadFloat ( &sigTemp ); 2191 else sigTemp = 0.0; 2192 } else { 2193 sigX = 0.0; 2194 sigY = 0.0; 2195 sigZ = 0.0; 2196 sigOcc = 0.0; 2197 sigTemp = 0.0; 2198 } 2199 2200 if (WhatIsSet & ASET_Anis_tFac) { 2201 f.ReadFloat ( &u11 ); 2202 f.ReadFloat ( &u22 ); 2203 f.ReadFloat ( &u33 ); 2204 f.ReadFloat ( &u12 ); 2205 f.ReadFloat ( &u13 ); 2206 f.ReadFloat ( &u23 ); 2207 if (WhatIsSet & ASET_Anis_tFSigma) { 2208 f.ReadFloat ( &su11 ); 2209 f.ReadFloat ( &su22 ); 2210 f.ReadFloat ( &su33 ); 2211 f.ReadFloat ( &su12 ); 2212 f.ReadFloat ( &su13 ); 2213 f.ReadFloat ( &su23 ); 2214 } else { 2215 su11 = 0.0; 2216 su22 = 0.0; 2217 su33 = 0.0; 2218 su12 = 0.0; 2219 su13 = 0.0; 2220 su23 = 0.0; 2221 } 2222 } else { 2223 u11 = 0.0; 2224 u22 = 0.0; 2225 u33 = 0.0; 2226 u12 = 0.0; 2227 u13 = 0.0; 2228 u23 = 0.0; 2229 su11 = 0.0; 2230 su22 = 0.0; 2231 su33 = 0.0; 2232 su12 = 0.0; 2233 su13 = 0.0; 2234 su23 = 0.0; 2235 } 2236 2237 f.ReadByte ( &nb ); 2238 if (nb>0) { 2239 Bond = new AtomBond[nb]; 2240 for (i=0;i<nb;i++) { 2241 f.ReadInt ( &k ); 2242 if (k>0) f.ReadByte ( &(Bond[i].order) ); 2243 else Bond[i].order = 0; 2244 // we place *index* of bonded atom temporary on the place 2245 // of its pointer, and the pointer will be calculated 2246 // after Residue::read calls _setBonds(..). 2247 memcpy ( &(Bond[i].atom),&k,4 ); 2248 } 2249 } 2250 nBonds = nb; 2251 nBonds = nBonds | (nBonds << 8); 2252 2253 } 2254 _setBonds(PPAtom A)2255 void Atom::_setBonds ( PPAtom A ) { 2256 int i,k,nb; 2257 nb = nBonds & 0x000000FF; 2258 for (i=0;i<nb;i++) { 2259 memcpy ( &k,&(Bond[i].atom),4 ); 2260 if (k>0) Bond[i].atom = A[k]; 2261 else Bond[i].atom = NULL; 2262 } 2263 } 2264 2265 MakeFactoryFunctions(Atom)2266 MakeFactoryFunctions(Atom) 2267 2268 2269 2270 // =========================== Residue =========================== 2271 2272 2273 void AtomStat::Init() { 2274 2275 nAtoms = 0; 2276 2277 xmin = MaxReal; xmax = MinReal; xm = 0.0; xm2 = 0.0; 2278 ymin = MaxReal; ymax = MinReal; ym = 0.0; ym2 = 0.0; 2279 zmin = MaxReal; zmax = MinReal; zm = 0.0; zm2 = 0.0; 2280 2281 occ_min = MaxReal; occ_max = MinReal; occ_m = 0.0; occ_m2 = 0.0; 2282 tFmin = MaxReal; tFmax = MinReal; tFm = 0.0; tFm2 = 0.0; 2283 2284 u11_min = MaxReal; u11_max = MinReal; u11_m = 0.0; u11_m2 = 0.0; 2285 u22_min = MaxReal; u22_max = MinReal; u22_m = 0.0; u22_m2 = 0.0; 2286 u33_min = MaxReal; u33_max = MinReal; u33_m = 0.0; u33_m2 = 0.0; 2287 u12_min = MaxReal; u12_max = MinReal; u12_m = 0.0; u12_m2 = 0.0; 2288 u13_min = MaxReal; u13_max = MinReal; u13_m = 0.0; u13_m2 = 0.0; 2289 u23_min = MaxReal; u23_max = MinReal; u23_m = 0.0; u23_m2 = 0.0; 2290 2291 WhatIsSet = ASET_All; 2292 2293 finished = false; 2294 2295 } 2296 Finish()2297 void AtomStat::Finish() { 2298 realtype v; 2299 2300 if (!finished) { 2301 2302 finished = true; 2303 2304 if (nAtoms>0) { 2305 2306 v = nAtoms; 2307 2308 xm /= v; xm2 /= v; 2309 ym /= v; ym2 /= v; 2310 zm /= v; zm2 /= v; 2311 2312 occ_m /= v; occ_m2 /= v; 2313 tFm /= v; tFm2 /= v; 2314 2315 u11_m /= v; u11_m2 /= v; 2316 u22_m /= v; u22_m2 /= v; 2317 u33_m /= v; u33_m2 /= v; 2318 u12_m /= v; u12_m2 /= v; 2319 u13_m /= v; u13_m2 /= v; 2320 u23_m /= v; u23_m2 /= v; 2321 } 2322 } 2323 2324 } 2325 GetMaxSize()2326 realtype AtomStat::GetMaxSize() { 2327 realtype r; 2328 r = RMax(xmax-xmin,ymax-ymin); 2329 r = RMax(r,zmax-zmin); 2330 return RMax(r,0.0); 2331 } 2332 2333 2334 // ---------------------------------------------------------------- 2335 2336 Residue()2337 Residue::Residue() : UDData() { 2338 InitResidue(); 2339 } 2340 Residue(PChain Chain_Owner)2341 Residue::Residue ( PChain Chain_Owner ) : UDData() { 2342 InitResidue(); 2343 if (Chain_Owner) 2344 Chain_Owner->AddResidue ( this ); 2345 } 2346 Residue(PChain Chain_Owner,const ResName resName,int sqNum,const InsCode ins)2347 Residue::Residue ( PChain Chain_Owner, 2348 const ResName resName, 2349 int sqNum, 2350 const InsCode ins ) : UDData() { 2351 InitResidue(); 2352 seqNum = sqNum; 2353 strcpy_css ( name,pstr(resName) ); 2354 strcpy_css ( insCode,pstr(ins) ); 2355 if (Chain_Owner) 2356 Chain_Owner->AddResidue ( this ); 2357 } 2358 Residue(io::RPStream Object)2359 Residue::Residue ( io::RPStream Object ) : UDData(Object) { 2360 InitResidue(); 2361 } 2362 ~Residue()2363 Residue::~Residue() { 2364 FreeMemory(); 2365 if (chain) chain->_ExcludeResidue ( name,seqNum,insCode ); 2366 } 2367 2368 InitResidue()2369 void Residue::InitResidue() { 2370 strcpy ( name ,"---" ); // residue name 2371 strcpy ( label_comp_id,"---" ); // assigned residue name 2372 label_asym_id[0] = char(0); // assigned chain Id 2373 seqNum = -MaxInt; // residue sequence number 2374 label_seq_id = -MaxInt; // assigned residue sequence number 2375 label_entity_id = 1; // assigned entity id 2376 strcpy ( insCode,"" ); // residue insertion code 2377 chain = NULL; // reference to chain 2378 index = -1; // undefined index in chain 2379 nAtoms = 0; // number of atoms in the residue 2380 AtmLen = 0; // length of atom array 2381 atom = NULL; // array of atoms 2382 Exclude = true; 2383 SSE = SSE_None; 2384 } 2385 SetChain(PChain Chain_Owner)2386 void Residue::SetChain ( PChain Chain_Owner ) { 2387 chain = Chain_Owner; 2388 } 2389 2390 GetResidueNo()2391 int Residue::GetResidueNo() { 2392 if (chain) return chain->GetResidueNo ( seqNum,insCode ); 2393 else return -1; 2394 } 2395 SetChainID(const ChainID chID)2396 void Residue::SetChainID ( const ChainID chID ) { 2397 if (chain) 2398 chain->SetChainID ( chID ); 2399 } 2400 2401 GetCenter(realtype & x,realtype & y,realtype & z)2402 int Residue::GetCenter ( realtype & x, realtype & y, 2403 realtype & z ) { 2404 int i,k; 2405 x = 0.0; 2406 y = 0.0; 2407 z = 0.0; 2408 k = 0; 2409 for (i=0;i<nAtoms;i++) 2410 if (atom[i]) { 2411 if (!atom[i]->Ter) { 2412 x += atom[i]->x; 2413 y += atom[i]->y; 2414 z += atom[i]->z; 2415 k++; 2416 } 2417 } 2418 if (k>0) { 2419 x /= k; 2420 y /= k; 2421 z /= k; 2422 return 0; 2423 } 2424 return 1; 2425 } 2426 GetCoordHierarchy()2427 void * Residue::GetCoordHierarchy() { 2428 if (chain) return chain->GetCoordHierarchy(); 2429 return NULL; 2430 } 2431 GetAltLocations(int & nAltLocs,PAltLoc & aLoc,rvector & occupancy,int & alflag)2432 void Residue::GetAltLocations ( int & nAltLocs, 2433 PAltLoc & aLoc, 2434 rvector & occupancy, 2435 int & alflag ) { 2436 int i,j,k, nal,nal1; 2437 realtype occ1; 2438 bool B; 2439 PAltLoc aL; 2440 rvector occ; 2441 bvector alv; 2442 2443 aLoc = NULL; 2444 occupancy = NULL; 2445 nAltLocs = 0; 2446 alflag = ALF_NoAltCodes; 2447 2448 if (nAtoms>0) { 2449 2450 // temporary array for altcodes 2451 aL = new AltLoc[nAtoms]; 2452 // temporary array for occupancies 2453 GetVectorMemory ( occ,nAtoms,0 ); 2454 // temporary array for checking altcodes 2455 GetVectorMemory ( alv,nAtoms,0 ); 2456 for (i=0;i<nAtoms;i++) 2457 alv[i] = false; 2458 2459 k = 0; // counts unique alternation codes 2460 nal = 0; 2461 for (i=0;i<nAtoms;i++) 2462 if (atom[i]) { 2463 if (!atom[i]->Ter) { 2464 // Find if the alternation code of ith atom is 2465 // a new one. 2466 B = false; 2467 for (j=0;(j<k) && (!B);j++) 2468 B = !strcmp(atom[i]->altLoc,aL[j]); 2469 if (!B) { 2470 // that's a new altcode, get its occupancy 2471 if (atom[i]->WhatIsSet & ASET_Occupancy) 2472 occ[k] = atom[i]->occupancy; 2473 else occ[k] = -1.0; 2474 // store new altcode in temporary array 2475 strcpy ( aL[k],atom[i]->altLoc ); 2476 // check consistency of the altcode data if: 2477 // a) the data was not found wrong so far 2478 // b) this atom name has not been checked before 2479 // c) altcode is not the "empty"-altcode 2480 if ((!(alflag & ALF_Mess)) && (!alv[i]) && 2481 (atom[i]->altLoc[0])) { 2482 B = false; // will be set true if "empty"-altcode 2483 // is found for current atom name 2484 nal1 = 0; // counts the number of different altcodes 2485 // for current atom name 2486 occ1 = 0.0; // will count the sum of occupancies for 2487 // current atom name 2488 for (j=0;j<nAtoms;j++) 2489 if (atom[j]) { 2490 if ((!atom[j]->Ter) && 2491 (!strcmp(atom[j]->name,atom[i]->name))) { 2492 if (atom[j]->WhatIsSet & ASET_Occupancy) 2493 occ1 += atom[j]->occupancy; 2494 if (!atom[j]->altLoc[0]) B = true; 2495 alv[j] = true; // mark it as "checked" 2496 nal1++; 2497 } 2498 } 2499 if (!(alflag & (ALF_EmptyAltLoc | ALF_NoEmptyAltLoc))) { 2500 if (B) alflag |= ALF_EmptyAltLoc; 2501 else alflag |= ALF_NoEmptyAltLoc; 2502 } else if (((alflag & ALF_EmptyAltLoc) && (!B)) || 2503 ((alflag & ALF_NoEmptyAltLoc) && (B))) 2504 alflag |= ALF_Mess; 2505 if ((occ[k]>=0) && (fabs(1.0-occ1)>0.01)) 2506 alflag |= ALF_Occupancy; 2507 if (nal==0) // first time just remember the number 2508 nal = nal1; // of different altcodes 2509 else if (nal!=nal1) // check if number of different altcodes 2510 alflag |= ALF_Mess; // is not the same through the residue 2511 } 2512 k++; 2513 } 2514 } 2515 } 2516 if (k>0) { 2517 aLoc = new AltLoc[k]; 2518 GetVectorMemory ( occupancy,k,0 ); 2519 for (i=0;i<k;i++) { 2520 strcpy ( aLoc[i],aL[i] ); 2521 occupancy[i] = occ[i]; 2522 } 2523 nAltLocs = k; 2524 } 2525 2526 delete[] aL; 2527 FreeVectorMemory ( occ,0 ); 2528 FreeVectorMemory ( alv,0 ); 2529 2530 } 2531 2532 } 2533 GetNofAltLocations()2534 int Residue::GetNofAltLocations() { 2535 int i,j,k; 2536 bool B; 2537 k = 0; 2538 for (i=0;i<nAtoms;i++) 2539 if (atom[i]) { 2540 if (!atom[i]->Ter) { 2541 B = false; 2542 for (j=0;(j<i) && (!B);j++) 2543 if (atom[j]) { 2544 if (!atom[j]->Ter) 2545 B = !strcmp(atom[i]->altLoc,atom[j]->altLoc); 2546 } 2547 if (!B) k++; 2548 } 2549 } 2550 return k; 2551 } 2552 SetResID(const ResName resName,int sqNum,const InsCode ins)2553 void Residue::SetResID ( const ResName resName, int sqNum, 2554 const InsCode ins ) { 2555 strcpy_css ( name,pstr(resName) ); 2556 seqNum = sqNum; 2557 strcpy_css ( insCode,pstr(ins) ); 2558 strcpy (label_comp_id,name ); 2559 } 2560 FreeMemory()2561 void Residue::FreeMemory() { 2562 // NOTE: individual atoms are disposed here as well! 2563 DeleteAllAtoms(); 2564 if (atom) delete[] atom; 2565 atom = NULL; 2566 nAtoms = 0; 2567 AtmLen = 0; 2568 } 2569 ExpandAtomArray(int nAdd)2570 void Residue::ExpandAtomArray ( int nAdd ) { 2571 int i; 2572 PPAtom atom1; 2573 AtmLen += abs(nAdd); 2574 atom1 = new PAtom[AtmLen]; 2575 for (i=0;i<nAtoms;i++) 2576 atom1[i] = atom[i]; 2577 for (i=nAtoms;i<AtmLen;i++) 2578 atom1[i] = NULL; 2579 if (atom) delete[] atom; 2580 atom = atom1; 2581 } 2582 _AddAtom(PAtom atm)2583 int Residue::_AddAtom ( PAtom atm ) { 2584 // Adds atom to the residue 2585 int i; 2586 for (i=0;i<nAtoms;i++) 2587 if (atom[i]==atm) return -i; // this atom is already there 2588 if (nAtoms>=AtmLen) 2589 ExpandAtomArray ( nAtoms+10-AtmLen ); 2590 atom[nAtoms] = atm; 2591 atom[nAtoms]->residue = this; 2592 nAtoms++; 2593 return 0; 2594 } 2595 AddAtom(PAtom atm)2596 int Residue::AddAtom ( PAtom atm ) { 2597 // AddAtom(..) adds atom to the residue. If residue is associated 2598 // with a coordinate hierarchy, and atom 'atm' is not, the latter 2599 // is checked in automatically. If atom 'atm' belongs to any 2600 // coordinate hierarchy (even though that of the residue), it is 2601 // *copied* rather than simply taken over, and is checked in. 2602 // If residue is not associated with a coordinate hierarchy, all 2603 // added atoms will be checked in automatically once the residue 2604 // is checked in. 2605 PRoot manager; 2606 PResidue res; 2607 int i; 2608 2609 for (i=0;i<nAtoms;i++) 2610 if (atom[i]==atm) return -i; // this atom is already there 2611 2612 if (nAtoms>=AtmLen) 2613 ExpandAtomArray ( nAtoms+10-AtmLen ); 2614 2615 if (atm->GetCoordHierarchy()) { 2616 atom[nAtoms] = newAtom(); 2617 atom[nAtoms]->Copy ( atm ); 2618 } else { 2619 res = atm->GetResidue(); 2620 if (res) 2621 for (i=0;i<res->nAtoms;i++) 2622 if (res->atom[i]==atm) { 2623 res->atom[i] = NULL; 2624 break; 2625 } 2626 atom[nAtoms] = atm; 2627 } 2628 2629 atom[nAtoms]->residue = this; 2630 manager = PRoot(GetCoordHierarchy()); 2631 if (manager) 2632 manager->CheckInAtom ( 0,atom[nAtoms] ); 2633 2634 nAtoms++; 2635 2636 return nAtoms; 2637 2638 } 2639 InsertAtom(PAtom atm,int position)2640 int Residue::InsertAtom ( PAtom atm, int position ) { 2641 // InsertAtom(..) inserts atom into the specified position of 2642 // the residue. If residue is associated with a coordinate hierarchy, 2643 // and atom 'atm' is not, the latter is checked in automatically. 2644 // If atom 'atm' belongs to any coordinate hierarchy (even though 2645 // that of the residue), it is *copied* rather than simply taken 2646 // over, and is checked in. 2647 // If residue is not associated with a coordinate hierarchy, all 2648 // added atoms will be checked in automatically once the residue 2649 // is checked in. 2650 PRoot manager; 2651 PResidue res; 2652 int i,pos; 2653 2654 for (i=0;i<nAtoms;i++) 2655 if (atom[i]==atm) return -i; // this atom is already there 2656 2657 if (nAtoms>=AtmLen) 2658 ExpandAtomArray ( nAtoms+10-AtmLen ); 2659 2660 pos = IMin(position,nAtoms); 2661 for (i=nAtoms;i>pos;i--) 2662 atom[i] = atom[i-1]; 2663 2664 if (atm->GetCoordHierarchy()) { 2665 atom[pos] = newAtom(); 2666 atom[pos]->Copy ( atm ); 2667 } else { 2668 res = atm->GetResidue(); 2669 if (res) 2670 for (i=0;i<res->nAtoms;i++) 2671 if (res->atom[i]==atm) { 2672 res->atom[i] = NULL; 2673 break; 2674 } 2675 atom[pos] = atm; 2676 } 2677 2678 atom[pos]->residue = this; 2679 manager = PRoot(GetCoordHierarchy()); 2680 if (manager) 2681 manager->CheckInAtom ( 0,atom[pos] ); 2682 2683 nAtoms++; 2684 2685 return nAtoms; 2686 2687 } 2688 InsertAtom(PAtom atm,const AtomName aname)2689 int Residue::InsertAtom ( PAtom atm, const AtomName aname ) { 2690 // This version inserts before the atom with given name. If such 2691 // name is not found, the atom is appended to the end. 2692 int i; 2693 i = 0; 2694 while (i<nAtoms) 2695 if (!atom[i]) i++; 2696 else if (!strcmp(aname,atom[i]->name)) break; 2697 else i++; 2698 return InsertAtom ( atm,i ); 2699 } 2700 2701 CheckInAtoms()2702 void Residue::CheckInAtoms() { 2703 PRoot manager; 2704 int i; 2705 manager = PRoot(GetCoordHierarchy()); 2706 if (manager) 2707 for (i=0;i<nAtoms;i++) 2708 if (atom[i]) { 2709 if (atom[i]->index<0) 2710 manager->CheckInAtom ( 0,atom[i] ); 2711 } 2712 } 2713 2714 _ExcludeAtom(int kndex)2715 int Residue::_ExcludeAtom ( int kndex ) { 2716 // deletes atom from the residue 2717 int i,k; 2718 2719 if (!Exclude) return 0; 2720 2721 k = -1; 2722 for (i=0;(i<nAtoms) && (k<0);i++) 2723 if (atom[i]) { 2724 if (atom[i]->index==kndex) k = i; 2725 } 2726 2727 if (k>=0) { 2728 for (i=k+1;i<nAtoms;i++) 2729 atom[i-1] = atom[i]; 2730 nAtoms--; 2731 } 2732 2733 if (nAtoms<=0) return 1; 2734 else return 0; 2735 2736 } 2737 2738 PDBASCIIAtomDump(io::RFile f)2739 void Residue::PDBASCIIAtomDump ( io::RFile f ) { 2740 int i; 2741 for (i=0;i<nAtoms;i++) 2742 if (atom[i]) 2743 atom[i]->PDBASCIIDump ( f ); 2744 } 2745 MakeAtomCIF(mmcif::PData CIF)2746 void Residue::MakeAtomCIF ( mmcif::PData CIF ) { 2747 int i; 2748 for (i=0;i<nAtoms;i++) 2749 if (atom[i]) 2750 atom[i]->MakeCIF ( CIF ); 2751 } 2752 2753 Copy(PResidue res)2754 void Residue::Copy ( PResidue res ) { 2755 // 2756 // Modify Residue::Copy and both Residues::_copy methods 2757 // simultaneously! 2758 // 2759 // This function will nake a copy of residue res in 'this' one. 2760 // All atoms are copied, none is moved regardless to the association 2761 // with coordinate hierarchy. If 'this' residue is associated with 2762 // a coordinate hierarchy, all atoms are checked in. 2763 PRoot manager; 2764 int i; 2765 2766 FreeMemory(); 2767 2768 seqNum = res->seqNum; 2769 label_seq_id = res->label_seq_id; 2770 label_entity_id = res->label_entity_id; 2771 index = res->index; 2772 AtmLen = res->nAtoms; 2773 SSE = res->SSE; 2774 strcpy ( name ,res->name ); 2775 strcpy ( label_comp_id,res->label_comp_id ); 2776 strcpy ( label_asym_id,res->label_asym_id ); 2777 strcpy ( insCode ,res->insCode ); 2778 2779 if (AtmLen>0) { 2780 atom = new PAtom[AtmLen]; 2781 nAtoms = 0; 2782 for (i=0;i<res->nAtoms;i++) 2783 if (res->atom[i]) { 2784 atom[nAtoms] = newAtom(); 2785 atom[nAtoms]->Copy ( res->atom[i] ); 2786 atom[nAtoms]->SetResidue ( this ); 2787 nAtoms++; 2788 } 2789 for (i=nAtoms;i<AtmLen;i++) 2790 atom[i] = NULL; 2791 manager = PRoot(GetCoordHierarchy()); 2792 if (manager) 2793 manager->CheckInAtoms ( 0,atom,nAtoms ); 2794 } 2795 2796 } 2797 2798 _copy(PResidue res)2799 void Residue::_copy ( PResidue res ) { 2800 // Modify both Residue::_copy and Residue::Copy methods 2801 // simultaneously! 2802 // 2803 // will work properly only if atomic arrays 2804 // this->chain->model->GetAtom() and 2805 // res->chain->model->GetAtom() are identical 2806 // 2807 int i; 2808 PPAtom A; 2809 2810 FreeMemory(); 2811 2812 seqNum = res->seqNum; 2813 label_seq_id = res->label_seq_id; 2814 label_entity_id = res->label_entity_id; 2815 index = res->index; 2816 nAtoms = res->nAtoms; 2817 SSE = res->SSE; 2818 strcpy ( name ,res->name ); 2819 strcpy ( label_comp_id,res->label_comp_id ); 2820 strcpy ( label_asym_id,res->label_asym_id ); 2821 strcpy ( insCode ,res->insCode ); 2822 2823 AtmLen = nAtoms; 2824 A = NULL; 2825 if (chain) { 2826 if (chain->model) 2827 A = chain->model->GetAllAtoms(); 2828 } 2829 if ((nAtoms>0) && (A)) { 2830 atom = new PAtom[nAtoms]; 2831 for (i=0;i<nAtoms;i++) { 2832 atom[i] = A[res->atom[i]->index-1]; 2833 atom[i]->SetResidue ( this ); 2834 } 2835 } else { 2836 nAtoms = 0; 2837 AtmLen = 0; 2838 } 2839 2840 } 2841 _copy(PResidue res,PPAtom atm,int & atom_index)2842 void Residue::_copy ( PResidue res, PPAtom atm, 2843 int & atom_index ) { 2844 // modify both Residue::_copy and Residue::Copy methods 2845 // simultaneously! 2846 // 2847 // This function physically copies the atoms, creating new atom 2848 // instances and putting them into array 'atm' sequentially from 2849 // 'atom_index' position. 'atom_index' is modified (advanced). 2850 // 2851 int i; 2852 2853 FreeMemory(); 2854 2855 seqNum = res->seqNum; 2856 label_seq_id = res->label_seq_id; 2857 label_entity_id = res->label_entity_id; 2858 index = res->index; 2859 nAtoms = res->nAtoms; 2860 SSE = res->SSE; 2861 strcpy ( name ,res->name ); 2862 strcpy ( label_comp_id,res->label_comp_id ); 2863 strcpy ( label_asym_id,res->label_asym_id ); 2864 strcpy ( insCode ,res->insCode ); 2865 2866 AtmLen = nAtoms; 2867 if (AtmLen>0) { 2868 atom = new PAtom[AtmLen]; 2869 for (i=0;i<nAtoms;i++) 2870 if (res->atom[i]) { 2871 if (!atm[atom_index]) atm[atom_index] = newAtom(); 2872 atm[atom_index]->Copy ( res->atom[i] ); 2873 atm[atom_index]->residue = this; 2874 atm[atom_index]->index = atom_index+1; 2875 atom[i] = atm[atom_index]; 2876 atom_index++; 2877 } else 2878 atom[i] = NULL; 2879 } 2880 2881 } 2882 2883 GetAtomStatistics(RAtomStat AS)2884 void Residue::GetAtomStatistics ( RAtomStat AS ) { 2885 AS.Init(); 2886 CalAtomStatistics ( AS ); 2887 AS.Finish(); 2888 } 2889 CalAtomStatistics(RAtomStat AS)2890 void Residue::CalAtomStatistics ( RAtomStat AS ) { 2891 // AS must be initialized. The function only accumulates 2892 // the statistics. 2893 int i; 2894 for (i=0;i<nAtoms;i++) 2895 if (atom[i]) 2896 atom[i]->CalAtomStatistics ( AS ); 2897 } 2898 2899 GetChain()2900 PChain Residue::GetChain() { 2901 return chain; 2902 } 2903 GetModel()2904 PModel Residue::GetModel() { 2905 if (chain) return (PModel)chain->model; 2906 else return NULL; 2907 } 2908 2909 GetModelNum()2910 int Residue::GetModelNum() { 2911 if (chain) { 2912 if (chain->model) 2913 return chain->model->GetSerNum(); 2914 } 2915 return 0; 2916 } 2917 GetChainID()2918 pstr Residue::GetChainID() { 2919 if (chain) return chain->chainID; 2920 return pstr(""); 2921 } 2922 GetLabelAsymID()2923 pstr Residue::GetLabelAsymID() { 2924 return label_asym_id; 2925 } 2926 GetResName()2927 pstr Residue::GetResName() { 2928 return name; 2929 } 2930 GetLabelCompID()2931 pstr Residue::GetLabelCompID() { 2932 return label_comp_id; 2933 } 2934 GetAASimilarity(const ResName resName)2935 int Residue::GetAASimilarity ( const ResName resName ) { 2936 return mmdb::GetAASimilarity ( pstr(name),pstr(resName) ); 2937 } 2938 GetAASimilarity(PResidue res)2939 int Residue::GetAASimilarity ( PResidue res ) { 2940 return mmdb::GetAASimilarity ( name,res->name ); 2941 } 2942 GetAAHydropathy()2943 realtype Residue::GetAAHydropathy() { 2944 return mmdb::GetAAHydropathy ( name ); 2945 } 2946 SetResName(const ResName resName)2947 void Residue::SetResName ( const ResName resName ) { 2948 strcpy ( name,resName ); 2949 } 2950 GetSeqNum()2951 int Residue::GetSeqNum() { 2952 return seqNum; 2953 } 2954 GetLabelSeqID()2955 int Residue::GetLabelSeqID() { 2956 return label_seq_id; 2957 } 2958 GetLabelEntityID()2959 int Residue::GetLabelEntityID() { 2960 return label_entity_id; 2961 } 2962 GetInsCode()2963 pstr Residue::GetInsCode() { 2964 return insCode; 2965 } 2966 isAminoacid()2967 bool Residue::isAminoacid () { 2968 return mmdb::isAminoacid ( name ); 2969 } 2970 isNucleotide()2971 bool Residue::isNucleotide() { 2972 return mmdb::isNucleotide ( name ); 2973 } 2974 isDNARNA()2975 int Residue::isDNARNA() { 2976 return mmdb::isDNARNA ( name ); 2977 } 2978 isSugar()2979 bool Residue::isSugar() { 2980 return mmdb::isSugar ( name ); 2981 } 2982 isSolvent()2983 bool Residue::isSolvent() { 2984 return mmdb::isSolvent ( name ); 2985 } 2986 isModRes()2987 bool Residue::isModRes() { 2988 PChain chn; 2989 PModRes modRes; 2990 int nModRes,i; 2991 chn = GetChain(); 2992 if (chn) { 2993 nModRes = chn->GetNofModResidues(); 2994 for (i=0;i<nModRes;i++) { 2995 modRes = chn->GetModResidue ( i ); 2996 if (modRes) { 2997 if ((!strcmp(modRes->resName,name)) && 2998 (modRes->seqNum==seqNum) && 2999 (!strcmp(modRes->insCode,insCode))) 3000 return true; 3001 } 3002 } 3003 3004 } 3005 return false; 3006 } 3007 isInSelection(int selHnd)3008 bool Residue::isInSelection ( int selHnd ) { 3009 PRoot manager = (PRoot)GetCoordHierarchy(); 3010 PMask mask; 3011 if (manager) { 3012 mask = manager->GetSelMask ( selHnd ); 3013 if (mask) return CheckMask ( mask ); 3014 } 3015 return false; 3016 } 3017 3018 isNTerminus()3019 bool Residue::isNTerminus() { 3020 PPResidue Res; 3021 int i,j,nRes; 3022 if (chain) { 3023 chain->GetResidueTable ( Res,nRes ); 3024 i = 0; 3025 j = -1; 3026 while ((i<nRes) && (j<0)) { 3027 if (Res[i]) j = i; 3028 i++; 3029 } 3030 if (j>=0) 3031 return (Res[j]->index==index); 3032 } 3033 return false; 3034 } 3035 isCTerminus()3036 bool Residue::isCTerminus() { 3037 PPResidue Res; 3038 int i,j,nRes; 3039 if (chain) { 3040 chain->GetResidueTable ( Res,nRes ); 3041 i = nRes-1; 3042 j = -1; 3043 while ((i>=0) && (j<0)) { 3044 if (Res[i]) j = i; 3045 i--; 3046 } 3047 if (j>=0) 3048 return (Res[j]->index==index); 3049 } 3050 return false; 3051 } 3052 3053 GetResidueID(pstr ResidueID)3054 pstr Residue::GetResidueID ( pstr ResidueID ) { 3055 ResidueID[0] = char(0); 3056 if (chain) { 3057 if (chain->model) 3058 sprintf ( ResidueID,"/%i/",chain->model->GetSerNum() ); 3059 else strcpy ( ResidueID,"/-/" ); 3060 strcat ( ResidueID,chain->chainID ); 3061 } else 3062 strcpy ( ResidueID,"/-/-" ); 3063 ParamStr ( ResidueID,pstr("/"),seqNum ); 3064 strcat ( ResidueID,"(" ); 3065 strcat ( ResidueID,name ); 3066 strcat ( ResidueID,")" ); 3067 if (insCode[0]) { 3068 strcat ( ResidueID,"." ); 3069 strcat ( ResidueID,insCode ); 3070 } 3071 return ResidueID; 3072 } 3073 3074 CheckID(int * snum,const InsCode inscode,const ResName resname)3075 int Residue::CheckID ( int * snum, 3076 const InsCode inscode, 3077 const ResName resname ) { 3078 if (snum) { 3079 if (*snum!=seqNum) return 0; 3080 } 3081 if (inscode) { 3082 if ((inscode[0]!='*') && (strcmp(inscode,insCode))) return 0; 3083 } 3084 if (!resname) return 1; 3085 if ((resname[0]!='*') && (strcmp(resname,name))) return 0; 3086 return 1; 3087 } 3088 CheckIDS(cpstr CID)3089 int Residue::CheckIDS ( cpstr CID ) { 3090 ChainID chn; 3091 InsCode inscode; 3092 ResName resname; 3093 AtomName atm; 3094 Element elm; 3095 AltLoc aloc; 3096 pstr p1,p2; 3097 int mdl,sn,rc; 3098 3099 rc = ParseAtomPath ( CID,mdl,chn,sn,inscode,resname, 3100 atm,elm,aloc,NULL ); 3101 // rc = ParseResID ( CID,sn,inscode,resname ); 3102 3103 if (rc>=0) { 3104 p1 = NULL; 3105 p2 = NULL; 3106 if (inscode[0]!='*') p1 = inscode; 3107 if (resname[0]!='*') p2 = resname; 3108 if (!rc) return CheckID ( &sn ,p1,p2 ); 3109 else return CheckID ( NULL,p1,p2 ); 3110 } 3111 return 0; 3112 3113 } 3114 3115 3116 // -------------------- Extracting atoms ------------------------- 3117 GetNumberOfAtoms()3118 int Residue::GetNumberOfAtoms() { 3119 return nAtoms; 3120 } 3121 GetNumberOfAtoms(bool countTers)3122 int Residue::GetNumberOfAtoms ( bool countTers ) { 3123 int i,na; 3124 na = 0; 3125 for (i=0;i<nAtoms;i++) 3126 if (atom[i]) { 3127 if (countTers || (!atom[i]->Ter)) na++; 3128 } 3129 return na; 3130 } 3131 GetAtom(const AtomName aname,const Element elname,const AltLoc aloc)3132 PAtom Residue::GetAtom ( const AtomName aname, 3133 const Element elname, 3134 const AltLoc aloc ) { 3135 int i; 3136 for (i=0;i<nAtoms;i++) 3137 if (atom[i]) { 3138 if (atom[i]->CheckID(aname,elname,aloc)) 3139 return atom[i]; 3140 } 3141 return NULL; 3142 } 3143 GetAtom(int atomNo)3144 PAtom Residue::GetAtom ( int atomNo ) { 3145 if ((0<=atomNo) && (atomNo<nAtoms)) 3146 return atom[atomNo]; 3147 return NULL; 3148 } 3149 GetAtomTable(PPAtom & atomTable,int & NumberOfAtoms)3150 void Residue::GetAtomTable ( PPAtom & atomTable, int & NumberOfAtoms ) { 3151 atomTable = atom; 3152 NumberOfAtoms = nAtoms; 3153 } 3154 GetAtomTable1(PPAtom & atomTable,int & NumberOfAtoms)3155 void Residue::GetAtomTable1 ( PPAtom & atomTable, int & NumberOfAtoms ) { 3156 int i,j; 3157 if (atomTable) delete[] atomTable; 3158 if (nAtoms>0) { 3159 atomTable = new PAtom[nAtoms]; 3160 j = 0; 3161 for (i=0;i<nAtoms;i++) 3162 if (atom[i]) { 3163 if (!atom[i]->Ter) 3164 atomTable[j++] = atom[i]; 3165 } 3166 NumberOfAtoms = j; 3167 } else { 3168 atomTable = NULL; 3169 NumberOfAtoms = 0; 3170 } 3171 } 3172 TrimAtomTable()3173 void Residue::TrimAtomTable() { 3174 int i,j; 3175 j = 0; 3176 for (i=0;i<nAtoms;i++) 3177 if (atom[i]) { 3178 if (j<i) { 3179 atom[j] = atom[i]; 3180 atom[i] = NULL; 3181 } 3182 j++; 3183 } 3184 nAtoms = j; 3185 } 3186 3187 3188 // --------------------- Deleting atoms -------------------------- 3189 DeleteAtom(const AtomName aname,const Element elname,const AltLoc aloc)3190 int Residue::DeleteAtom ( const AtomName aname, 3191 const Element elname, 3192 const AltLoc aloc ) { 3193 // apply Root::FinishStructEdit() after all editings are done! 3194 // returns number of deleted atoms 3195 int i,k,nA,kndex; 3196 PPAtom A; 3197 3198 A = NULL; 3199 nA = 0; 3200 if (chain) { 3201 if (chain->model) { 3202 A = chain->model->GetAllAtoms(); 3203 nA = chain->model->GetNumberOfAllAtoms(); 3204 } 3205 } 3206 3207 k = 0; 3208 for (i=0;i<nAtoms;i++) 3209 if (atom[i]) { 3210 if (atom[i]->CheckID(aname,elname,aloc)) { 3211 k++; 3212 kndex = atom[i]->index; 3213 if ((0<kndex) && (kndex<=nA)) A[kndex-1] = NULL; 3214 Exclude = false; 3215 delete atom[i]; 3216 atom[i] = NULL; 3217 Exclude = true; 3218 } 3219 } 3220 3221 return k; 3222 3223 } 3224 DeleteAtom(int atomNo)3225 int Residue::DeleteAtom ( int atomNo ) { 3226 // apply Root::FinishStructEdit() after all editings are done! 3227 // returns number of deleted atoms 3228 int kndex,nA; 3229 PPAtom A; 3230 3231 if ((0<=atomNo) && (atomNo<nAtoms)) { 3232 if (atom[atomNo]) { 3233 A = NULL; 3234 nA = 0; 3235 if (chain) { 3236 if (chain->model) { 3237 A = chain->model->GetAllAtoms(); 3238 nA = chain->model->GetNumberOfAllAtoms(); 3239 } 3240 } 3241 kndex = atom[atomNo]->index; 3242 if ((0<kndex) && (kndex<=nA)) A[kndex-1] = NULL; 3243 Exclude = false; 3244 delete atom[atomNo]; 3245 atom[atomNo] = NULL; 3246 Exclude = true; 3247 return 1; 3248 } 3249 } 3250 3251 return 0; 3252 3253 } 3254 3255 DeleteAllAtoms()3256 int Residue::DeleteAllAtoms() { 3257 int i,k,nA,kndex; 3258 PPAtom A; 3259 3260 Exclude = false; 3261 3262 A = NULL; 3263 nA = 0; 3264 if (chain) { 3265 if (chain->model) { 3266 A = chain->model->GetAllAtoms(); 3267 nA = chain->model->GetNumberOfAllAtoms(); 3268 } 3269 } 3270 3271 k = 0; 3272 for (i=0;i<nAtoms;i++) 3273 if (atom[i]) { 3274 k++; 3275 kndex = atom[i]->index; 3276 if ((0<kndex) && (kndex<=nA)) A[kndex-1] = NULL; 3277 delete atom[i]; 3278 atom[i] = NULL; 3279 } 3280 nAtoms = 0; 3281 3282 Exclude = true; 3283 3284 return k; 3285 3286 } 3287 3288 DeleteAltLocs()3289 int Residue::DeleteAltLocs() { 3290 // This function leaves only alternative location with maximal 3291 // occupancy, if those are equal or unspecified, the one with 3292 // "least" alternative location indicator. 3293 // The function returns the number of deleted atoms. The atom 3294 // table remains untrimmed, so that nAtoms are wrong until that 3295 // is done. Tables are trimmed by FinishStructEdit() or 3296 // explicitely. 3297 PPAtom A; 3298 AtomName aname; 3299 AltLoc aLoc,aL; 3300 realtype occupancy,occ; 3301 int nA,i,i1,i2,j,k,n,kndex; 3302 3303 A = NULL; 3304 nA = 0; 3305 if (chain) { 3306 if (chain->model) { 3307 A = chain->model->GetAllAtoms(); 3308 nA = chain->model->GetNumberOfAllAtoms(); 3309 } 3310 } 3311 Exclude = false; 3312 3313 n = 0; 3314 for (i=0;i<nAtoms;i++) 3315 3316 if (atom[i]) { 3317 if (!atom[i]->Ter) { 3318 occupancy = atom[i]->GetOccupancy(); 3319 strcpy ( aname,atom[i]->name ); 3320 strcpy ( aLoc ,atom[i]->altLoc ); 3321 i1 = -1; 3322 i2 = i; 3323 k = 0; 3324 for (j=i+1;j<nAtoms;j++) 3325 if (atom[j]) { 3326 if ((!atom[j]->Ter) && (!strcmp(atom[j]->name,aname))) { 3327 k++; 3328 occ = atom[j]->GetOccupancy(); 3329 if (occ>occupancy) { 3330 occupancy = occ; 3331 i1 = j; 3332 } 3333 if (aLoc[0]) { 3334 strcpy ( aL,atom[j]->altLoc ); 3335 if (!aL[0]) { 3336 aLoc[0] = char(0); 3337 i2 = j; 3338 } else if (strcmp(aL,aLoc)<0) { 3339 strcpy ( aLoc,aL ); 3340 i2 = j; 3341 } 3342 } 3343 } 3344 } 3345 if (k>0) { 3346 if (i1<0) { 3347 if (atom[i]->WhatIsSet & ASET_Occupancy) i1 = i; 3348 else i1 = i2; 3349 } 3350 for (j=i;j<nAtoms;j++) 3351 if ((j!=i1) && atom[j]) { 3352 if ((!atom[j]->Ter) && (!strcmp(atom[j]->name,aname))) { 3353 n++; 3354 kndex = atom[j]->index; 3355 if ((0<kndex) && (kndex<=nA)) A[kndex-1] = NULL; 3356 delete atom[j]; 3357 atom[j] = NULL; 3358 } 3359 } 3360 } 3361 } 3362 } 3363 3364 Exclude = true; 3365 3366 return n; 3367 3368 } 3369 ApplyTransform(const mat44 & TMatrix)3370 void Residue::ApplyTransform ( const mat44 & TMatrix ) { 3371 // transforms all coordinates by multiplying with matrix TMatrix 3372 int i; 3373 for (i=0;i<nAtoms;i++) 3374 if (atom[i]) { 3375 if (!atom[i]->Ter) 3376 atom[i]->Transform ( TMatrix ); 3377 } 3378 } 3379 3380 3381 3382 // ----------------------------------------------------------------- 3383 3384 MaskAtoms(PMask Mask)3385 void Residue::MaskAtoms ( PMask Mask ) { 3386 int i; 3387 for (i=0;i<nAtoms;i++) 3388 if (atom[i]) atom[i]->SetMask ( Mask ); 3389 } 3390 UnmaskAtoms(PMask Mask)3391 void Residue::UnmaskAtoms ( PMask Mask ) { 3392 int i; 3393 for (i=0;i<nAtoms;i++) 3394 if (atom[i]) atom[i]->RemoveMask ( Mask ); 3395 } 3396 3397 3398 3399 // ------- user-defined data handlers 3400 PutUDData(int UDDhandle,int iudd)3401 int Residue::PutUDData ( int UDDhandle, int iudd ) { 3402 if (UDDhandle & UDRF_RESIDUE) 3403 return UDData::putUDData ( UDDhandle,iudd ); 3404 else return UDDATA_WrongUDRType; 3405 } 3406 PutUDData(int UDDhandle,realtype rudd)3407 int Residue::PutUDData ( int UDDhandle, realtype rudd ) { 3408 if (UDDhandle & UDRF_RESIDUE) 3409 return UDData::putUDData ( UDDhandle,rudd ); 3410 else return UDDATA_WrongUDRType; 3411 } 3412 PutUDData(int UDDhandle,cpstr sudd)3413 int Residue::PutUDData ( int UDDhandle, cpstr sudd ) { 3414 if (UDDhandle & UDRF_RESIDUE) 3415 return UDData::putUDData ( UDDhandle,sudd ); 3416 else return UDDATA_WrongUDRType; 3417 } 3418 GetUDData(int UDDhandle,int & iudd)3419 int Residue::GetUDData ( int UDDhandle, int & iudd ) { 3420 if (UDDhandle & UDRF_RESIDUE) 3421 return UDData::getUDData ( UDDhandle,iudd ); 3422 else return UDDATA_WrongUDRType; 3423 } 3424 GetUDData(int UDDhandle,realtype & rudd)3425 int Residue::GetUDData ( int UDDhandle, realtype & rudd ) { 3426 if (UDDhandle & UDRF_RESIDUE) 3427 return UDData::getUDData ( UDDhandle,rudd ); 3428 else return UDDATA_WrongUDRType; 3429 } 3430 GetUDData(int UDDhandle,pstr sudd,int maxLen)3431 int Residue::GetUDData ( int UDDhandle, pstr sudd, int maxLen ) { 3432 if (UDDhandle & UDRF_RESIDUE) 3433 return UDData::getUDData ( UDDhandle,sudd,maxLen ); 3434 else return UDDATA_WrongUDRType; 3435 } 3436 GetUDData(int UDDhandle,pstr & sudd)3437 int Residue::GetUDData ( int UDDhandle, pstr & sudd ) { 3438 if (UDDhandle & UDRF_RESIDUE) 3439 return UDData::getUDData ( UDDhandle,sudd ); 3440 else return UDDATA_WrongUDRType; 3441 } 3442 3443 3444 #define NOmaxdist2 12.25 3445 isMainchainHBond(PResidue res)3446 bool Residue::isMainchainHBond ( PResidue res ) { 3447 // Test if there is main chain Hbond between PCRes1 (donor) and 3448 // PCRes2 (acceptor). 3449 // As defined Kabsch & Sanders 3450 // This probably needs the option of supporting alternative criteria 3451 PAtom NAtom,OAtom,Atom; 3452 realtype abx,aby,abz; 3453 realtype acx,acy,acz; 3454 realtype bcx,bcy,bcz; 3455 realtype absq,acsq,bcsq; 3456 3457 NAtom = GetAtom ( "N" ); 3458 OAtom = res->GetAtom ( "O" ); 3459 Atom = res->GetAtom ( "C" ); 3460 3461 if (NAtom && OAtom && Atom) { 3462 3463 abx = OAtom->x - NAtom->x; 3464 aby = OAtom->y - NAtom->y; 3465 abz = OAtom->z - NAtom->z; 3466 absq = abx*abx + aby*aby + abz*abz; 3467 3468 3469 if (absq<=NOmaxdist2) { 3470 3471 acx = NAtom->x - Atom->x; 3472 acy = NAtom->y - Atom->y; 3473 acz = NAtom->z - Atom->z; 3474 3475 bcx = Atom->x - OAtom->x; 3476 bcy = Atom->y - OAtom->y; 3477 bcz = Atom->z - OAtom->z; 3478 3479 acsq = acx*acx + acy*acy + acz*acz; 3480 bcsq = bcx*bcx + bcy*bcy + bcz*bcz; 3481 3482 return (acos((bcsq+absq-acsq)/(2.0*sqrt(bcsq*absq)))>=Pi/2.0); 3483 3484 } 3485 3486 } 3487 3488 return false; 3489 3490 } 3491 3492 write(io::RFile f)3493 void Residue::write ( io::RFile f ) { 3494 int i; 3495 bool shortBinary = false; 3496 byte Version=3; 3497 3498 f.WriteByte ( &Version ); 3499 3500 if (nAtoms>0) { 3501 if (atom[0]->WhatIsSet & ASET_CompactBinary) 3502 shortBinary = true; 3503 } 3504 3505 f.WriteBool ( &shortBinary ); 3506 3507 if (!shortBinary) { 3508 UDData::write ( f ); 3509 f.WriteInt ( &label_seq_id ); 3510 f.WriteInt ( &label_entity_id ); 3511 f.WriteTerLine ( label_comp_id,false ); 3512 f.WriteTerLine ( label_asym_id,false ); 3513 } 3514 3515 f.WriteInt ( &seqNum ); 3516 f.WriteInt ( &index ); 3517 f.WriteInt ( &nAtoms ); 3518 f.WriteByte ( &SSE ); 3519 f.WriteTerLine ( name ,false ); 3520 f.WriteTerLine ( insCode,false ); 3521 for (i=0;i<nAtoms;i++) 3522 f.WriteInt ( &(atom[i]->index) ); 3523 3524 } 3525 read(io::RFile f)3526 void Residue::read ( io::RFile f ) { 3527 // IMPORTANT: array Atom in Root class should be 3528 // read prior calling this function! 3529 PPAtom A; 3530 int i,k; 3531 byte Version; 3532 bool shortBinary; 3533 3534 FreeMemory(); 3535 3536 f.ReadByte ( &Version ); 3537 f.ReadBool ( &shortBinary ); 3538 3539 if (!shortBinary) { 3540 UDData::read ( f ); 3541 f.ReadInt ( &label_seq_id ); 3542 f.ReadInt ( &label_entity_id ); 3543 f.ReadTerLine ( label_comp_id,false ); 3544 f.ReadTerLine ( label_asym_id,false ); 3545 } 3546 3547 f.ReadInt ( &seqNum ); 3548 f.ReadInt ( &index ); 3549 f.ReadInt ( &nAtoms ); 3550 f.ReadByte ( &SSE ); 3551 f.ReadTerLine ( name ,false ); 3552 f.ReadTerLine ( insCode,false ); 3553 AtmLen = nAtoms; 3554 A = NULL; 3555 if (chain) { 3556 if (chain->model) 3557 A = chain->model->GetAllAtoms(); 3558 } 3559 if ((nAtoms>0) && (A)) { 3560 atom = new PAtom[nAtoms]; 3561 for (i=0;i<nAtoms;i++) { 3562 f.ReadInt ( &k ); 3563 atom[i] = A[k-1]; 3564 atom[i]->SetResidue ( this ); 3565 atom[i]->_setBonds ( A ); 3566 } 3567 } else { 3568 for (i=0;i<nAtoms;i++) 3569 f.ReadInt ( &k ); 3570 nAtoms = 0; 3571 AtmLen = 0; 3572 } 3573 3574 } 3575 3576 3577 /* === old function, keep for a while 3578 void Residue::read ( io::RFile f ) { 3579 // IMPORTANT: array Atom in Root class should be 3580 // read prior calling this function! 3581 PPAtom A; 3582 int i,k; 3583 byte Version; 3584 3585 FreeMemory (); 3586 3587 UDData::read ( f ); 3588 3589 f.ReadByte ( &Version ); 3590 f.ReadInt ( &seqNum ); 3591 if (Version>1) { 3592 f.ReadInt ( &label_seq_id ); 3593 f.ReadInt ( &label_entity_id ); 3594 } 3595 f.ReadInt ( &index ); 3596 f.ReadInt ( &nAtoms ); 3597 f.ReadByte ( &SSE ); 3598 f.ReadTerLine ( name,false ); 3599 if (Version>1) { 3600 f.ReadTerLine ( label_comp_id,false ); 3601 f.ReadTerLine ( label_asym_id,false ); 3602 } 3603 f.ReadTerLine ( insCode,false ); 3604 AtmLen = nAtoms; 3605 A = NULL; 3606 if (chain) { 3607 if (chain->model) 3608 A = chain->model->GetAllAtoms(); 3609 } 3610 if ((nAtoms>0) && (A)) { 3611 atom = new PAtom[nAtoms]; 3612 for (i=0;i<nAtoms;i++) { 3613 f.ReadInt ( &k ); 3614 atom[i] = A[k-1]; 3615 atom[i]->SetResidue ( this ); 3616 atom[i]->_setBonds ( A ); 3617 } 3618 } else { 3619 for (i=0;i<nAtoms;i++) 3620 f.ReadInt ( &k ); 3621 nAtoms = 0; 3622 AtmLen = 0; 3623 } 3624 } 3625 */ 3626 3627 MakeFactoryFunctions(Residue) 3628 3629 } // namespace mmdb 3630