00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #ifndef __MSD_CHEM_H__
00021 #include "msd_chem.h"
00022 #endif
00023
00024 #ifdef __WITH_XMMDB__
00025
00026 HRESULT MSDTitle::LoadWarehouse(msdbConnect& db, int nDepId)
00027 { char select_stm[3072];
00028 int nDepDepId = nDepId;
00029 otl_datetime otl_item;
00030 struct SPdbTitle
00031 {
00032 char m_szAccessionCode[10];
00033 char m_szExpTitle[256];
00034 char m_szTitle[2001];
00035 char m_szDate[30];
00036
00037 } title;
00038 struct SPdbtitleInd
00039 {
00040 short m_nAccessionCode;
00041 short m_nExpTitle;
00042 short m_nTitle;
00043 short m_nDate;
00044
00045 } ind;
00046 sprintf(select_stm,"select /*+ INDEX(ENTRY)*/ ACCESSION_CODE, PDB_EXPERIMENT_TYPE, TITLE, CREATION_DATE "
00047 "from ENTRY where ENTRY_ID = :nDepDepId <int>");
00048 try{
00049 msdbStream imsd (50, select_stm, db);
00050 imsd << nDepDepId;
00051
00052 while (!imsd.eof()){
00053 imsd >> title.m_szAccessionCode;
00054
00055 ind.m_nAccessionCode =imsd.is_null();
00056 imsd >> title.m_szExpTitle;
00057 ind.m_nExpTitle=imsd.is_null();
00058 imsd >> title.m_szTitle;
00059 ind.m_nTitle=imsd.is_null();
00060 imsd >> otl_item;
00061 ind.m_nDate=imsd.is_null();
00062
00063 if (!ind.m_nAccessionCode)
00064 {
00065 trimRight(title.m_szAccessionCode);
00066 strncpy(idCode, title.m_szAccessionCode, sizeof(idCode));
00067 idCode[sizeof(idCode) - 1] = 0;
00068 }
00069 if (!ind.m_nExpTitle)
00070 {
00071 trimRight(title.m_szExpTitle);
00072 ExpData.AddData(new CExpData(title.m_szExpTitle));
00073 }
00074 if (!ind.m_nTitle)
00075 {
00076 trimRight(title.m_szTitle);
00077 Title.AddData(new CTitleLine(title.m_szTitle));
00078 }
00079
00080 if (! ind.m_nDate)
00081 {sprintf(title.m_szDate ,"%d-%d-%d",otl_item.year,otl_item.month,otl_item.day);
00082 }
00083 }
00084 }
00085 catch(msdbException& p){
00086 cerr<<"select title error"<<endl;
00087 cerr<<p.msg<<endl;
00088 cerr<<p.stm_text<<endl;
00089 cerr<<p.var_info<<endl;
00090 }
00091 return S_OK;
00092 }
00093
00094 HRESULT MSDCryst::LoadWarehouse(msdbConnect& db, int nDepId)
00095 { char select_stm[3072];
00096 HRESULT hr = S_OK;
00097 int nDepDepId = nDepId;
00098 struct SCryst
00099 {
00100 double m_nId, m_dblAlpha, m_dblBeta, m_dblGamma, m_dblA, m_dblB, m_dblC;
00101
00102 } cr;
00103 struct SCrystInd
00104 {
00105 short m_nId, m_nAlpha, m_nBeta, m_nGamma, m_nA, m_nB, m_nC;
00106
00107 } ind;
00108
00109 sprintf(select_stm, "select /*+ INDEX(XRAY_DATA)*/ ENTRY_ID, ANGLE_ALPHA, ANGLE_BETA, ANGLE_GAMMA, "
00110 "LENGTH_A, LENGTH_B, LENGTH_C from XRAY_DATA where ENTRY_ID = :nDepDepId<int> ");
00111
00112
00113 try{
00114
00115 msdbStream imsd (50, select_stm, db);
00116 imsd << nDepDepId;
00117 while (!imsd.eof()){
00118 imsd >> cr.m_nId;
00119 ind.m_nId = imsd.is_null();
00120 imsd >> cr.m_dblAlpha;
00121 ind.m_nAlpha = imsd.is_null();
00122 imsd >> cr.m_dblBeta;
00123 ind.m_nBeta = imsd.is_null();
00124 imsd >> cr.m_dblGamma;
00125 ind.m_nGamma = imsd.is_null();
00126 imsd >> cr.m_dblA;
00127 ind.m_nA = imsd.is_null();
00128 imsd >> cr.m_dblB;
00129 ind.m_nB = imsd.is_null();
00130 imsd >> cr.m_dblC;
00131 ind.m_nC = imsd.is_null();
00132 if (!(ind.m_nAlpha || ind.m_nBeta || ind.m_nGamma))
00133 {
00134 alpha = cr.m_dblAlpha;
00135 beta = cr.m_dblBeta;
00136 gamma = cr.m_dblGamma;
00137 WhatIsSet |= CSET_CellParams1;
00138 }
00139 if (!(ind.m_nA || ind.m_nB || ind.m_nC))
00140 {
00141 a = cr.m_dblA;
00142 b = cr.m_dblB;
00143 c = cr.m_dblC;
00144 WhatIsSet |= CSET_CellParams2;
00145 }
00146 WhatIsSet |= CSET_SpaceGroup;
00147 this->CalcCoordTransforms();
00148 }
00149
00150 }
00151 catch(msdbException& p){
00152 cerr<<"get cryst info"<<endl;
00153 cerr<<p.msg<<endl;
00154 cerr<<p.stm_text<<endl;
00155 cerr<<p.var_info<<endl;
00156 }
00157 return hr;
00158 }
00159
00160 HRESULT MSDManager::LoadWarehouse(msdbConnect& db, int nDepId)
00161 {
00162 char select_stm[3072];
00163 int nDepDepId = nDepId;
00164 struct SChain
00165 {
00166
00167 int m_nAssemblyId;
00168 int m_nAssemblySerial;
00169
00170 int m_nGenAssemblyId;
00171 int m_nSymOpId;
00172 double m_T11, m_T12, m_T13, m_T14,
00173 m_T21, m_T22, m_T23, m_T24,
00174 m_T31, m_T32, m_T33, m_T34;
00175
00176 int m_nChainId;
00177 char m_szChainPdbId[2];
00178 char m_szChainType[4];
00179 char m_chMsdCode;
00180 } chain;
00181 struct SChainInd
00182 {
00183
00184 short m_nAssemblyId, m_nAssemblySerial;
00185
00186 short m_nGenAssemblyId, m_nSymOpId,
00187 m_T11, m_T12, m_T13, m_T14,
00188 m_T21, m_T22, m_T23, m_T24,
00189 m_T31, m_T32, m_T33, m_T34;
00190
00191 short m_nChainId, m_nChainPdbId, m_nChainType, m_nMsdCode;
00192 } chainInd;
00193 struct SData
00194 {
00195 int m_nGenAssemblyId;
00196
00197 int m_nModelId;
00198 int m_nModelSerial;
00199
00200 int m_nResidueId;
00201 int m_nResiduePdbSeq;
00202 char m_szResidueCode[12];
00203 char m_szResiduePdbInsCode[2];
00204 int m_nSupChemCompId;
00205
00206 int m_nAtomId;
00207 int m_nAtomSiteId;
00208 int m_nSupChemAtomId;
00209 int m_nAtomSerial;
00210 char m_szElement[8];
00211 int m_nCharge;
00212 char m_szAtomName[12];
00213 char m_szAltLoc[10];
00214 double m_dblAltId;
00215 double m_dblOccupancy;
00216 double m_dblDrX, m_dblDrY, m_dblDrZ;
00217 int m_nOrdering;
00218 int m_nAssemblyId;
00219 } data;
00220 struct SDataInd
00221 {
00222 short m_nGenAssemblyId;
00223
00224 short m_nModelId,
00225 m_nModelSerial;
00226
00227 short m_nResidueId,
00228 m_nResiduePdbSeq,
00229 m_nResidueCode,
00230 m_nResiduePdbInsCode,
00231 m_nSupChemCompId;
00232
00233 short m_nAtomId,
00234 m_nAtomSiteId,
00235 m_nSupChemAtomId,
00236 m_nAtomSerial,
00237 m_nElement,
00238 m_nCharge,
00239 m_nAtomName,
00240 m_nAltLoc,
00241 m_nAltId,
00242 m_nOccupancy,
00243 m_nDrX, m_nDrY, m_nDrZ,
00244 m_nOrdering,
00245 m_nAssemblyId;
00246 } dataInd;
00247 struct SPlaneAtom
00248 {
00249 int m_nSupChemPlaneId;
00250 int m_nSupChemCompId;
00251 int m_nSupChemAtomId;
00252 } planeAtom;
00253 struct SPlaneAtomInd
00254 {
00255 short
00256 m_nSupChemPlaneId,
00257 m_nSupChemCompId,
00258 m_nSupChemAtomId;
00259 } planeAtomInd;
00260 struct SPlaneCenter
00261 {
00262 int m_nId;
00263 int m_nPlaneId;
00264 int m_nComponentId;
00265 int m_nChainId;
00266 int m_nModelId;
00267 double m_dblX, m_dblY, m_dblZ, m_dblNX, m_dblNY, m_dblNZ;
00268 char m_szAltLoc[10];
00269 double m_dblOccupancy;
00270 } planeCenter;
00271 struct SPlaneCenterInd
00272 {
00273 short m_nId,
00274 m_nPlaneId,
00275 m_nComponentId,
00276 m_nChainId,
00277 m_nModelId,
00278 m_nX, m_nY, m_nZ, m_nNX, m_nNY, m_nNZ,
00279 m_nAltLoc, m_nOccupancy;
00280 } planeCenterInd;
00281
00282
00283
00284 ResetManager ();
00285 FreeFileMemory();
00286 FType = MMDB_FILE_PDB;
00287 ignoreSegID = (Flags & MMDBF_IgnoreSegID )!=0;
00288 ignoreElement = (Flags & MMDBF_IgnoreElement)!=0;
00289 ignoreCharge = (Flags & MMDBF_IgnoreCharge )!=0;
00290 ignoreNonCoorPDBErrors = (Flags & MMDBF_IgnoreNonCoorPDBErrors)!=0;
00291
00292
00293 if (FAILED(((MSDTitle*)&Title)->LoadWarehouse(db, nDepId)))
00294 return E_FAIL;
00295 if (FAILED(((MSDCryst*)&Cryst)->LoadWarehouse(db, nDepId)))
00296 return E_FAIL;
00297
00298
00299 sprintf(select_stm,"select /*+ INDEX (CHAIN)*/"
00300 "CHAIN.ASSEMBLY_ID, CHAIN.ASSEMBLY_SERIAL, "
00301 "CHAIN.CHAIN_ID, CHAIN.INT_RSE_ID, "
00302 "MATRIX.M11, MATRIX.M12, MATRIX.M13, MATRIX.M14, "
00303 "MATRIX.M21, MATRIX.M22, MATRIX.M23, MATRIX.M24, "
00304 "MATRIX.M31, MATRIX.M32, MATRIX.M33, MATRIX.M34, "
00305 "CHAIN.CHAIN_ID, CHAIN.PDB_CODE, CHAIN.CHAIN_TYPE, "
00306 "substr(CHAIN.MSD_CODE, 0, 1)"
00307 "from CHAIN, MATRIX "
00308 "where MATRIX.MATRIX_ID (+) = CHAIN.MATRIX_ID AND "
00309 "CHAIN.ENTRY_ID = :nDepDepId<int> and CHAIN.CHAIN_ID > 0 ");
00310
00311 std::vector<CDBAssembly*> vpAssembly;
00312 CDBListChain* pLChain = NULL;
00313 CDBListChainAssembly* pLChainAssembly = NULL;
00314
00315
00316
00317 try{
00318
00319 msdbStream imsd (50, select_stm, db);
00320 imsd << nDepDepId;
00321 while (!imsd.eof()){
00322
00323 imsd >> chain.m_nAssemblyId;
00324 chainInd.m_nAssemblyId =imsd.is_null();
00325 imsd >> chain.m_nAssemblySerial;
00326 chainInd.m_nAssemblySerial=imsd.is_null();
00327
00328 imsd >> chain.m_nGenAssemblyId;
00329 chainInd.m_nGenAssemblyId=imsd.is_null();
00330 imsd >> chain.m_nSymOpId;
00331 chainInd.m_nSymOpId=imsd.is_null();
00332
00333 imsd >> chain.m_T11;
00334 chainInd.m_T11=imsd.is_null();
00335 imsd >> chain.m_T12;
00336 chainInd.m_T12=imsd.is_null();
00337 imsd >> chain.m_T13;
00338 chainInd.m_T13=imsd.is_null();
00339 imsd >> chain.m_T14;
00340 chainInd.m_T14=imsd.is_null();
00341 imsd >> chain.m_T21;
00342 chainInd.m_T21=imsd.is_null();
00343 imsd >> chain.m_T22;
00344 chainInd.m_T22=imsd.is_null();
00345 imsd >> chain.m_T23;
00346 chainInd.m_T23=imsd.is_null();
00347 imsd >> chain.m_T24;
00348 chainInd.m_T24=imsd.is_null();
00349 imsd >> chain.m_T31;
00350 chainInd.m_T31=imsd.is_null();
00351 imsd >> chain.m_T32;
00352 chainInd.m_T32=imsd.is_null();
00353 imsd >> chain.m_T33;
00354 chainInd.m_T33=imsd.is_null();
00355 imsd >> chain.m_T34;
00356 chainInd.m_T34=imsd.is_null();
00357
00358 imsd >> chain.m_nChainId;
00359 chainInd.m_nChainId=imsd.is_null();
00360 imsd >> chain.m_szChainPdbId;
00361 chainInd.m_nChainPdbId=imsd.is_null();
00362 imsd >> chain.m_szChainType;
00363 chainInd.m_nChainType=imsd.is_null();
00364 imsd >> chain.m_chMsdCode;
00365 chainInd.m_nMsdCode=imsd.is_null();
00366
00367 if (chainInd.m_nGenAssemblyId)
00368 continue;
00369 CDBAssembly* pAssembly = NULL;
00370 for (std::vector<CDBAssembly*>::const_iterator it = vpAssembly.begin();
00371 !pAssembly && it!= vpAssembly.end(); it++)
00372 {
00373 int nAssemblyId = 0;
00374 if ((*it)->GetId(&nAssemblyId) == S_OK && nAssemblyId == chain.m_nAssemblyId)
00375 pAssembly = *it;
00376 }
00377 if (!pAssembly)
00378 {
00379 pAssembly = new CDBAssembly(this, chain.m_nAssemblyId, chain.m_nAssemblySerial);
00380 vpAssembly.push_back(pAssembly);
00381 }
00382 CDBListChainAssembly* pChainAssembly = pLChainAssembly;
00383 for(; pChainAssembly; pChainAssembly = pChainAssembly->m_pPrev)
00384 {
00385 int nId;
00386 if (pChainAssembly->GetId(&nId) == S_OK && nId == chain.m_nGenAssemblyId)
00387 break;
00388 }
00389 if (!pChainAssembly)
00390 {
00391 pChainAssembly = new CDBListChainAssembly(pAssembly, chain.m_nGenAssemblyId, chain.m_nSymOpId, chain.m_nChainId,
00392 chain.m_T11, chain.m_T12, chain.m_T13, chain.m_T14,
00393 chain.m_T21, chain.m_T22, chain.m_T23, chain.m_T24,
00394 chain.m_T31, chain.m_T32, chain.m_T33, chain.m_T34,
00395 pLChainAssembly);
00396 pLChainAssembly = pChainAssembly;
00397 pAssembly->IncChainAssemblyCount();
00398 }
00399 CDBListChain* pChain = pLChain;
00400 for(;pChain; pChain = pChain->m_pPrev)
00401 {
00402 int nId;
00403 if (pChain->GetId(&nId) == S_OK && nId == chain.m_nChainId)
00404 break;
00405 }
00406 if (!pChain)
00407 {
00408 if (chainInd.m_nChainId)
00409 continue;
00410 if (!chainInd.m_nChainPdbId)
00411 trimRight(chain.m_szChainPdbId);
00412 else
00413 chain.m_szChainPdbId[0] = 0;
00414 pChain = new CDBListChain(chain.m_szChainPdbId, chain.m_nChainId, chain.m_szChainType, chain.m_chMsdCode, pLChain);
00415 pLChain = pChain;
00416 }
00417
00418 }
00419
00420 }
00421 catch(msdbException& p){
00422 cerr<<"get chain info"<<endl;
00423 cerr<<p.msg<<endl;
00424 cerr<<p.stm_text<<endl;
00425 cerr<<p.var_info<<endl;
00426 }
00427
00428
00429
00430
00431
00432
00433 sprintf(select_stm, "select /*+ INDEX(ATOM_DATA ATD_ENTRY_IND) */ "
00434 " CHAIN_ID, "
00435 " MODEL_ID, MODEL_SERIAL, "
00436 " RESIDUE_ID, "
00437 " RESIDUE_PDB_SEQ, "
00438 " nvl(CODE_3_LETTER, comp_code.id(CHEM_COMP_CODE)), "
00439 " RESIDUE_PDB_INSERT_CODE, CHEM_COMP_ID, "
00440 " ATOM_DATA_ID, ATOM_ID, CHEM_ATOM_ID, "
00441 " SERIAL, ELEMENT_SYMBOL, PDB_CHARGE, "
00442 " CHEM_ATOM_NAME, ALT_ID, "
00443 " OCCUPANCY, "
00444 " X, Y, Z, CHEM_ATOM_ORDERING, ASSEMBLY_ID "
00445 " from ATOM_DATA "
00446 " where ENTRY_ID = :nDepDepId <int>");
00447
00448
00449
00450 std::vector<CDBModel*> vpModel;
00451 CDBListResidue* pLResidue = NULL;
00452 CDBListAtom* pLAtom = NULL;
00453 CDBListModelAssembly* pLModelAssembly = NULL;
00454 CDBListChain* pLModelChain = NULL;
00455 std::vector<int> vnSupChemCompId;
00456 vnSupChemCompId.reserve(500);
00457 m_vpSupChemCompGraph.reserve(500);
00458
00459 try{
00460 msdbStream imsd (50, select_stm, db);
00461 imsd << nDepDepId;
00462 while (!imsd.eof()){
00463
00464 imsd >> data.m_nGenAssemblyId;
00465 dataInd.m_nGenAssemblyId=imsd.is_null();
00466
00467 imsd >> data.m_nModelId;
00468 dataInd.m_nModelId=imsd.is_null();
00469 imsd >> data.m_nModelSerial;
00470 dataInd.m_nModelSerial=imsd.is_null();
00471
00472 imsd >> data.m_nResidueId;
00473 dataInd.m_nResidueId=imsd.is_null();
00474 imsd >> data.m_nResiduePdbSeq;
00475 dataInd.m_nResiduePdbSeq=imsd.is_null();
00476 imsd >> data.m_szResidueCode;
00477 dataInd.m_nResidueCode=imsd.is_null();
00478 imsd >> data.m_szResiduePdbInsCode;
00479 dataInd.m_nResiduePdbInsCode=imsd.is_null();
00480 imsd >> data.m_nSupChemCompId;
00481 dataInd.m_nSupChemCompId=imsd.is_null();
00482
00483 imsd >> data.m_nAtomId;
00484 dataInd.m_nAtomId=imsd.is_null();
00485 imsd >> data.m_nAtomSiteId;
00486 dataInd.m_nAtomSiteId=imsd.is_null();
00487 imsd >> data.m_nSupChemAtomId;
00488 dataInd.m_nSupChemAtomId=imsd.is_null();
00489 imsd >> data.m_nAtomSerial;
00490 dataInd.m_nAtomSerial=imsd.is_null();
00491 imsd >> data.m_szElement;
00492 dataInd.m_nElement=imsd.is_null();
00493 imsd >> data.m_nCharge;
00494 dataInd.m_nCharge=imsd.is_null();
00495 imsd >> data.m_szAtomName;
00496 dataInd.m_nAtomName=imsd.is_null();
00497 imsd >> data.m_dblAltId;
00498 dataInd.m_nAltId=imsd.is_null();
00499 imsd >> data.m_dblOccupancy;
00500 dataInd.m_nOccupancy=imsd.is_null();
00501 imsd >> data.m_dblDrX;
00502 dataInd.m_nDrX=imsd.is_null();
00503 imsd >> data.m_dblDrY;
00504 dataInd.m_nDrY=imsd.is_null();
00505 imsd >> data.m_dblDrZ;
00506 dataInd.m_nDrZ=imsd.is_null();
00507 imsd >> data.m_nOrdering;
00508 dataInd.m_nOrdering=imsd.is_null();
00509 imsd >> data.m_nAssemblyId;
00510 dataInd.m_nAssemblyId=imsd.is_null();
00512
00513 if (dataInd.m_nAtomId) {
00514 cout << "wrong atom_id indicator" << endl;
00515 continue;
00516 }
00517 if (data.m_nAssemblyId < 0)
00518 continue;
00519
00520 CDBModel* pModel = NULL;
00521 for (std::vector<CDBModel*>::const_iterator it = vpModel.begin();
00522 !pModel && it != vpModel.end(); it++)
00523 {
00524 int nModelId = 0;
00525 if ((*it)->GetId(&nModelId) == S_OK && nModelId == data.m_nModelId)
00526 pModel = *it;
00527 }
00528 if (!pModel)
00529 {
00530 pModel = new CDBModel(this, data.m_nModelSerial, data.m_nModelId);
00531 vpModel.push_back(pModel);
00532 }
00533
00534 CDBListChainAssembly* pChainAssembly = pLChainAssembly;
00535 for(; pChainAssembly; pChainAssembly = pChainAssembly->m_pPrev)
00536 {
00537 int nId;
00538 if (pChainAssembly->GetId(&nId) == S_OK && nId == data.m_nGenAssemblyId)
00539 break;
00540 }
00541 CDBAssembly* pAssembly = NULL;
00542 if (!pChainAssembly) {
00543 cout << "absent chain:" << data.m_nGenAssemblyId << endl;
00544 continue;
00545 }
00546 pChainAssembly->GetAssembly(&pAssembly);
00547 int nAssemblyId = 0;
00548 pAssembly->GetId(&nAssemblyId);
00549 CDBListChain* pChain = pLModelChain;
00550 for(;pChain; pChain = pChain->m_pPrev)
00551 {
00552 int nId;
00553 if (pChain->GetModel() == pModel &&
00554 pChain->GetId(&nId) == S_OK && nId == data.m_nGenAssemblyId)
00555 break;
00556 }
00557 if (!pChain)
00558 {
00559 CDBListChain* pExChain = pLChain;
00560 int nId = 0;
00561 for(;pExChain; pExChain = pExChain->m_pPrev)
00562 {
00563 if (pExChain->GetId(&nId) == S_OK && nId == data.m_nGenAssemblyId)
00564 break;
00565 }
00566 char szType[40];
00567 pExChain->GetType(szType);
00568 pChain = new CDBListChain(pExChain->GetChainID(), nId, szType, pExChain->IsProtein() == S_OK ? 'P' : 'N', pLModelChain);
00569 pLModelChain = pChain;
00570 pChain->SetModel(pModel);
00571 pModel->IncChains();
00572 }
00573
00574 CDBListModelAssembly* pModelAssembly = pLModelAssembly;
00575 for(; pModelAssembly; pModelAssembly = pModelAssembly->m_pPrev)
00576 {
00577 CDBModel* pM = NULL;
00578 CDBAssembly* pA = NULL;
00579 if (pModelAssembly->GetModel(&pM) == S_OK &&
00580 pModelAssembly->GetAssembly(&pA) == S_OK &&
00581 pA == pAssembly && pM == pModel)
00582 break;
00583 }
00584 if (!pModelAssembly)
00585 {
00586 pModelAssembly = new CDBListModelAssembly(pModel, pAssembly, this, pLModelAssembly);
00587 pLModelAssembly = pModelAssembly;
00588 pAssembly->IncModelAssemblyCount();
00589 }
00590
00591 CDBListResidue* pResidue = pLResidue;
00592 for(;pResidue; pResidue = pResidue->m_pPrev)
00593 {
00594 int nId;
00595 if (pResidue->GetChain() == pChain &&
00596 pResidue->GetId(&nId) == S_OK && nId == data.m_nResidueId)
00597 break;
00598 }
00599 if (!pResidue)
00600 {
00601 char szType[40];
00602 pChain->GetType(szType);
00603 bool bIsHet = !(szType[0] == 'C' || szType[1] == 'C');
00604 if (!dataInd.m_nResidueCode)
00605 {
00606 trimRight(data.m_szResidueCode);
00607 if (islower(data.m_szResidueCode[0]))
00608 data.m_szResidueCode[0] = 0;
00609 else if (strlen(data.m_szResidueCode) > 3)
00610 data.m_szResidueCode[1] = 0;
00611 }
00612 else
00613 data.m_szResidueCode[0] = 0;
00614 if (!dataInd.m_nResiduePdbInsCode)
00615 trimRight(data.m_szResiduePdbInsCode);
00616 else
00617 data.m_szResiduePdbInsCode[0] = 0;
00618 if (dataInd.m_nResiduePdbSeq)
00619 data.m_nResiduePdbSeq = 0;
00620 if (dataInd.m_nSupChemCompId)
00621 data.m_nSupChemCompId = 0;
00622 pResidue = new CDBListResidue(data.m_szResidueCode, data.m_nResiduePdbSeq,
00623 data.m_szResiduePdbInsCode, data.m_nResidueId, bIsHet,
00624 data.m_nSupChemCompId, nAssemblyId, pLResidue);
00625 pLResidue = pResidue;
00626 pResidue->SetChain(pChain);
00627 pResidue->index = ++m_nResidueIndex;
00628 pResidue->nAtoms = 0;
00629 pChain->IncResidues();
00630 if (data.m_nSupChemCompId)
00631 {
00632 std::vector<int>::iterator it = std::find(vnSupChemCompId.begin(), vnSupChemCompId.end(), data.m_nSupChemCompId);
00633 if (it == vnSupChemCompId.end())
00634 vnSupChemCompId.push_back(data.m_nSupChemCompId);
00635 }
00636 }
00637 if (!dataInd.m_nElement)
00638 trimRight(data.m_szElement);
00639 if (!dataInd.m_nAtomName)
00640 trimRight(data.m_szAtomName);
00641 if (!dataInd.m_nAltLoc)
00642 trimRight(data.m_szAltLoc);
00643 char szEmpty[2] = {0};
00644 char szCharge[10] = {0};
00645 if (!dataInd.m_nCharge)
00646 sprintf(szCharge, "%i3", data.m_nCharge);
00647 if (dataInd.m_nSupChemAtomId)
00648 data.m_nSupChemAtomId = 0;
00649 pLAtom = new CDBListAtom(data.m_nAtomId, data.m_nAtomSiteId, data.m_nSupChemAtomId,
00650 pLAtom);
00651
00652 pLAtom->SetAtomName(m_nAtomIndex + 1,
00653 dataInd.m_nOrdering ? m_nAtomIndex + 1 : data.m_nOrdering,
00654 dataInd.m_nAtomName ? szEmpty : data.m_szAtomName,
00655 dataInd.m_nAltLoc ? szEmpty : data.m_szAltLoc,
00656 szEmpty,
00657 dataInd.m_nElement ? szEmpty : data.m_szElement,
00658 szCharge);
00659 pLAtom->SetCoordinates(dataInd.m_nDrX ? 0 : data.m_dblDrX,
00660 dataInd.m_nDrY ? 0 : data.m_dblDrY,
00661 dataInd.m_nDrZ ? 0 : data.m_dblDrZ,
00662 dataInd.m_nOccupancy ? 0 : data.m_dblOccupancy,
00663 0);
00664 bool bIsHet = false;
00665 pResidue->IsHet(&bIsHet);
00666 pLAtom->Het = bIsHet;
00667 pLAtom->SetResidue(pResidue);
00668
00669 pResidue->nAtoms++;
00670 m_nAtomIndex++;
00671 }
00672 }
00673 catch(msdbException& p){
00674 cerr<<"get data info"<<endl;
00675 cerr<<p.msg<<endl;
00676 cerr<<p.stm_text<<endl;
00677 cerr<<p.var_info<<endl;
00678 }
00679
00680 sprintf(select_stm,"select /*+ INDEX_COMBINE(P PLANE_CENTER_ENTRY_ID PLANE_CENTER_ALT_ID) INDEX(ALT ALT_PK) */"
00681 " P.ID, P.PLANE_ID, P.RESIDUE_ID, P.CHAIN_ID, P.MODEL_ID, "
00682 " P.X, P.Y, P.Z, P.NX, P.NY, P.NZ, ALT.ALT_CODE, P.OCCUPANCY "
00683 " from PLANE_CENTER P, ALT "
00684 " where P.ENTRY_ID = :nDepDepId <int> and ALT.ALT_ID (+) = P.ALT_ID");
00685
00686
00687 int nPlaneCenterCount = 0;
00688
00689 try{
00690
00691 msdbStream imsd (50, select_stm, db);
00692 imsd << nDepDepId;
00693 while (!imsd.eof()){
00694
00695 imsd >> planeCenter.m_nId;
00696 planeCenterInd.m_nId=imsd.is_null();
00697 imsd >> planeCenter.m_nPlaneId;
00698 planeCenterInd.m_nPlaneId=imsd.is_null();
00699 imsd >> planeCenter.m_nComponentId;
00700 planeCenterInd.m_nComponentId=imsd.is_null();
00701 imsd >> planeCenter.m_nChainId;
00702 planeCenterInd.m_nChainId=imsd.is_null();
00703 imsd >> planeCenter.m_nModelId;
00704 planeCenterInd.m_nModelId=imsd.is_null();
00705 imsd >> planeCenter.m_dblX;
00706 planeCenterInd.m_nX=imsd.is_null();
00707 imsd >> planeCenter.m_dblY;
00708 planeCenterInd.m_nY=imsd.is_null();
00709 imsd >> planeCenter.m_dblZ;
00710 planeCenterInd.m_nZ=imsd.is_null();
00711 imsd >> planeCenter.m_dblNX;
00712 planeCenterInd.m_nNX=imsd.is_null();
00713 imsd >> planeCenter.m_dblNY;
00714 planeCenterInd.m_nNY=imsd.is_null();
00715 imsd >> planeCenter.m_dblNZ;
00716 planeCenterInd.m_nNZ=imsd.is_null();
00717 imsd >> planeCenter.m_szAltLoc;
00718 planeCenterInd.m_nAltLoc=imsd.is_null();
00719 imsd >> planeCenter.m_dblOccupancy;
00720 planeCenterInd.m_nOccupancy=imsd.is_null();
00721
00722 if (planeCenterInd.m_nId) {
00723 cout << "wrong plane indicator:" << planeCenterInd.m_nId << endl;
00724 continue;
00725 }
00726
00727 int nModelId = 0;
00728 int nChainId = 0;
00729 int nResidueId = 0;
00730 CDBListResidue *pResidue = pLResidue;
00731 for(; pResidue; pResidue = pResidue->m_pPrev)
00732 if (((CDBModel*)pResidue->GetChain()->GetModel())->GetId(&nModelId) == S_OK &&
00733 nModelId == planeCenter.m_nModelId &&
00734 ((CDBChain*)pResidue->GetChain())->GetId(&nChainId) == S_OK &&
00735 nChainId == planeCenter.m_nChainId &&
00736 pResidue->GetId(&nResidueId) == S_OK &&
00737 nResidueId == planeCenter.m_nComponentId)
00738 break;
00739 if (!pResidue) {
00740 cout << "residue is missed componentId:" << planeCenter.m_nComponentId << endl
00741 << "modelId:" << planeCenter.m_nModelId << endl
00742 << "chainId:" << planeCenter.m_nChainId << endl;
00743 continue;
00744 }
00745 if (planeCenterInd.m_nAltLoc)
00746 planeCenter.m_szAltLoc[0] = 0;
00747 else
00748 trimRight(planeCenter.m_szAltLoc);
00749 algebra::vec3 center, norm;
00750 SegID segId = {0};
00751 center[0] = planeCenter.m_dblX;
00752 center[1] = planeCenter.m_dblY;
00753 center[2] = planeCenter.m_dblZ;
00754 norm[0] = planeCenter.m_dblNX;
00755 norm[1] = planeCenter.m_dblNY;
00756 norm[2] = planeCenter.m_dblNZ;
00757 pLAtom = new CDBListPlaneAtom(planeCenter.m_nId, norm, pLAtom);
00758 char szEmpty[2] = {0};
00759 char szCharge[10] = "0";
00760 pLAtom->SetAtomName(m_nAtomIndex + 1,
00761 m_nAtomIndex + 1,
00762 "GVC",
00763 planeCenter.m_szAltLoc,
00764 segId,
00765 szEmpty,
00766 szCharge);
00767 pLAtom->SetCoordinates(center[0],
00768 center[1],
00769 center[2],
00770 planeCenter.m_dblOccupancy,
00771 0);
00772 bool bIsHet = false;
00773 pResidue->IsHet(&bIsHet);
00774 pLAtom->Het = bIsHet;
00775 pLAtom->SetResidue(pResidue);
00776 pResidue->nAtoms++;
00777 m_nAtomIndex++;
00778
00779 }
00780 }
00781 catch(msdbException& p){
00782 cerr<<"get planeCenter info"<<endl;
00783 cerr<<p.msg<<endl;
00784 cerr<<p.stm_text<<endl;
00785 cerr<<p.var_info<<endl;
00786 }
00787
00788
00789
00790 for(CDBListResidue* pResidue = pLResidue; pResidue; pResidue = pResidue->m_pPrev)
00791 {
00792 int nCount = pResidue->nAtoms;
00793 CDBListAtom* pAtom = pLAtom;
00794 pResidue->nAtoms = 0;
00795 pResidue->atom = new CAtom*[nCount];
00796 for(; pAtom; pAtom = pAtom->m_pPrev)
00797 {
00798 if (pAtom->GetResidue() == pResidue)
00799 pResidue->atom[pResidue->nAtoms++] = pAtom;
00800 }
00801 if (pResidue->nAtoms > 1)
00802 std::stable_sort(pResidue->atom, pResidue->atom + pResidue->nAtoms, atom_comp);
00803 }
00804
00805 m_nAssemblyCount = 0;
00806 m_vpAssembly = new CDBAssembly*[vpAssembly.size()];
00807 for (std::vector<CDBAssembly*>::iterator it = vpAssembly.begin();
00808 it != vpAssembly.end();
00809 it++)
00810 {
00811 m_vpAssembly[m_nAssemblyCount++] = *it;
00812 int nChainAssemblyCount = 0;
00813 int nModelAssemblyCount = 0;
00814 (*it)->GetChainAssemblyCount(&nChainAssemblyCount);
00815 (*it)->GetModelAssemblyCount(&nModelAssemblyCount);
00816 (*it)->PopulateChainAssemblies(pLChainAssembly, nChainAssemblyCount);
00817 (*it)->PopulateModelAssemblies(pLModelAssembly, nModelAssemblyCount);
00818 }
00819 for(CDBListChain* pChain = pLModelChain; pChain; pChain = pChain->m_pPrev)
00820 pChain->PopulateResidues(pLResidue, pChain->GetNumberOfResidues());
00821
00822 Model = new CModel*[vpModel.size()];
00823 nModels = 0;
00824 for (std::vector<CDBModel*>::iterator it = vpModel.begin();
00825 it != vpModel.end();
00826 it++)
00827 {
00828 Model[nModels++] = *it;
00829 (*it)->PopulateChains(pLModelChain, (*it)->GetNumberOfChains());
00830 }
00831
00832
00833 AtmLen = 0;
00834 nAtoms = m_nAtomIndex;
00835 Atom = new CAtom*[nAtoms];
00836 for (int i = 0; i < nModels; i++)
00837 {
00838 PCModel pModel = Model[i];
00839 for (int j = 0; j < pModel->GetNumberOfChains(); j++)
00840 {
00841 CDBChain* pChain = (CDBChain*)pModel->GetChain(j);
00842 char szType[12];
00843 pChain->GetType(szType);
00844
00845
00846 for (int k = 0; k < pChain->GetNumberOfResidues(); k++)
00847 {
00848 PCResidue pRes = pChain->GetResidue(k);
00849 for (int l = 0; l < pRes->nAtoms; l++)
00850 Atom[AtmLen++] = pRes->atom[l];
00851 }
00852 }
00853 }
00854 for (int i = 0; i < nAtoms; i++)
00855 ((CAtomWrapper*)Atom[i])->SetIndex(i + 1);
00856
00857
00858 for (std::vector<int>::const_iterator it = vnSupChemCompId.begin();
00859 it != vnSupChemCompId.end(); it++)
00860 m_vpSupChemCompGraph.push_back(buildGraphFromWarehouseComponent(db,*it, false));
00861
00862 for(CDBListChain* pExChain = pLChain; pExChain;)
00863 {
00864 CDBListChain* pChain = pExChain;
00865 pExChain = pExChain->m_pPrev;
00866 delete pChain;
00867 }
00868 PDBCleanup ( PDBCLEAN_ATNAME | PDBCLEAN_CHAIN_ORDER );
00869 return S_OK;
00870 }
00871
00872
00873
00875 PCCalcGraph buildGraphFromWarehouseComponent(msdbConnect& db, int nscc_id, bool bEliminateH)
00876 {
00877
00878 char select_stm[3072];
00879 int nSccId;
00880 struct SAtom
00881 {
00882 int m_nId;
00883 char m_szName[13];
00884 char m_szElement[16];
00885 double m_dblCharge;
00886 } atom;
00887 struct SAtomInd
00888 {
00889 short m_nId,m_nName,m_nElement, m_nCharge;
00890 } atomInd;
00891 struct SAtomPos
00892 {
00893 int m_nId;
00894 char m_szType[256];
00895 double m_dblX, m_dblY, m_dblZ;
00896 } atomPos;
00897 struct SAtomPosInd
00898 {
00899 short m_nId, m_nType, m_nX, m_nY, m_nZ;
00900 } atomPosInd;
00901 struct SBond
00902 {
00903 char m_szOrder[5];
00904 int m_nId1, m_nId2;
00905 } bond;
00906 struct SBondInd
00907 {
00908 short m_nOrder,m_nId1, m_nId2;
00909 } bondInd;
00910 char szEliminateType[3];
00911
00912 nSccId = nscc_id;
00913 strcpy(szEliminateType,bEliminateH ? "H" : " ");
00914 PCCalcGraph pGraph = new CCalcGraph(nSccId);
00915
00916
00917
00918
00919
00920 sprintf(select_stm, "select /*+ RULE */ a.chem_atom_id, a.name, e.symbol, a.charge "
00921 "from chem_atom a, element e "
00922 "where :nSccId <int> in (a.chem_comp_id) and "
00923 "a.element_id = e.element_id and e.symbol != :szEliminateType <char[3]> ");
00924
00925 try{
00926
00927 msdbStream imsd (50, select_stm, db);
00928 imsd << nSccId;
00929 imsd << szEliminateType;
00930 int nVertexId = 1;
00931 while (!imsd.eof()){
00932
00933
00934 imsd >> atom.m_nId;
00935 atomInd.m_nId=imsd.is_null();
00936 imsd >> atom.m_szName;
00937 atomInd.m_nName=imsd.is_null();
00938 imsd >> atom.m_szElement;
00939 atomInd.m_nElement=imsd.is_null();
00940 imsd >> atom.m_dblCharge;
00941 atomInd.m_nCharge=imsd.is_null();
00942 if (atomInd.m_nId)
00943 continue;
00944 trimRight(atom.m_szName);
00945 trimRight(atom.m_szElement);
00946 if (atomInd.m_nCharge)
00947 atom.m_dblCharge = 0;
00948 PCCalcVertex pVer = new CCalcVertex(pGraph, atom.m_szElement, atom.m_szName, nVertexId, atom.m_nId);
00949 pVer->SetCharge(atom.m_dblCharge);
00950 pGraph->AddVertex(pVer);
00951 nVertexId++;
00952 }
00953 }
00954 catch(msdbException& p){
00955 cerr<<"get atom info error"<<endl;
00956 cerr<<p.msg<<endl;
00957 cerr<<p.stm_text<<endl;
00958 cerr<<p.var_info<<endl;
00959 }
00960
00961
00962
00963 sprintf(select_stm,"select /*+ RULE */ a.chem_atom_id, rl.name, l.model_x, l.model_y, l.model_z "
00964 " from chem_atom a, chem_model_coord l, data_source rl, element e "
00965 " where :nSccId <int> in (a.chem_comp_id) and "
00966 " a.element_id = e.element_id and e.symbol != :szEliminateType <char[3]> and "
00967 " l.chem_atom_id = a.chem_atom_id and l.chem_model_coord_id = rl.data_source_id");
00968
00969
00970
00971
00972 try{
00973
00974 msdbStream imsd (50, select_stm, db);
00975 imsd << nSccId;
00976 imsd << szEliminateType;
00977
00978 while (!imsd.eof()){
00979 imsd >> atomPos.m_nId;
00980 atomPosInd.m_nId=imsd.is_null();
00981 imsd >> atomPos.m_szType;
00982 atomPosInd.m_nType=imsd.is_null();
00983 imsd >> atomPos.m_dblX;
00984 atomPosInd.m_nX=imsd.is_null();
00985 imsd >> atomPos.m_dblY;
00986 atomPosInd.m_nY=imsd.is_null();
00987 imsd >> atomPos.m_dblZ;
00988 atomPosInd.m_nZ=imsd.is_null();
00989
00990
00991 if (atomPosInd.m_nId)
00992 continue;
00993
00994 trimRight(atomPos.m_szType);
00995 PCCalcVertex pVertex = pGraph->GetVertexByUserId(atomPos.m_nId);
00996 if (!pVertex)
00997 throw "Vertex is null";
00998 pVertex->SetCoords(atomPos.m_szType, atomPos.m_dblX, atomPos.m_dblY, atomPos.m_dblZ);
00999
01000 }
01001 }
01002
01003 catch(msdbException& p){
01004 cerr<<"get atom position error"<<endl;
01005 cerr<<p.msg<<endl;
01006 cerr<<p.stm_text<<endl;
01007 cerr<<p.var_info<<endl;
01008 }
01009 catch (char * str)
01010 {
01011 cerr<<str<<endl;
01012 }
01013
01014
01015
01016
01017
01018
01019 sprintf(select_stm, "select /*+ RULE */ b.chem_bond_type, b.chem_atom_1_id, b.chem_atom_2_id "
01020 "from chem_bond b, chem_atom a1, "
01021 " chem_atom a2, element e1, element e2 "
01022 " where a1.chem_atom_id = b.chem_atom_1_id and a2.chem_atom_id = b.chem_atom_2_id and "
01023 " :nSccId <int> in (a1.chem_comp_id) and "
01024 " :nSccId in (a2.chem_comp_id) and "
01025 " a1.element_id = e1.element_id and a2.element_id = e2.element_id and "
01026 " e1.symbol != :szEliminateType <char[3]> and e2.symbol != :szEliminateType");
01027 try{
01028
01029 msdbStream imsd (50, select_stm, db);
01030 imsd << nSccId;
01031 imsd << szEliminateType;
01032
01033 while (!imsd.eof()){
01034
01035
01036 imsd >> bond.m_szOrder;
01037 bondInd.m_nOrder=imsd.is_null();
01038 imsd >> bond.m_nId1;
01039 bondInd.m_nId1=imsd.is_null();
01040 imsd >> bond.m_nId2;
01041 bondInd.m_nId2=imsd.is_null();
01042
01043
01044 if (bondInd.m_nId1 || bondInd.m_nId2)
01045 continue;
01046 trimRight(bond.m_szOrder);
01047 int nBondOrder = BOND_SINGLE;
01048 if (!strncmp(bond.m_szOrder, "sing", 4))
01049 nBondOrder = BOND_SINGLE;
01050 else if (!strncmp(bond.m_szOrder, "doub", 4))
01051 nBondOrder = BOND_DOUBLE;
01052 else if (!strncmp(bond.m_szOrder, "arom", 4))
01053 nBondOrder = BOND_AROMATIC;
01054 else if (!strncmp(bond.m_szOrder, "trip", 4))
01055 nBondOrder = BOND_TRIPLE;
01056 int nVertex1 = pGraph->GetVertexByUserId(bond.m_nId1)->GetID();
01057 int nVertex2 = pGraph->GetVertexByUserId(bond.m_nId2)->GetID();
01058 pGraph->AddEdge(new CCalcEdge(nVertex1, nVertex2, nBondOrder));
01059
01060 }
01061 }
01062 catch(msdbException& p){
01063 cerr<<"get bond error"<<endl;
01064 cerr<<p.msg<<endl;
01065 cerr<<p.stm_text<<endl;
01066 cerr<<p.var_info<<endl;
01067 }
01068
01069
01070 pGraph->Build(True);
01071 return pGraph;
01072 }
01074
01075 MakeFactoryFunctions(CAtomWrapper)
01076 MakeFactoryFunctions(CDBAtom)
01077 HRESULT msdGetEntry(msdbConnect& db, int nDepId, int nMaxCells, double dblMaxDist,
01078 connectivity_d* connect_d, int nDCount,
01079 connectivity_a* connect_a, int nACount)
01080 {
01081
01082 char select_stm[3072];
01083 int nDepDepId = nDepId;
01084 char szExpType[255] = {0};
01085 short nExpTypeInd = 0;
01086 short nInd;
01087 int nContactType;
01088 int nContactType1;
01089 int nContactType2;
01090 int nValueType;
01091 double dblValue, dblAngle, dblDist1, dblDist2;
01092 int nActiveSiteType1, nActiveSiteType2;
01093 int nAtomId1, nAtomId2;
01094 char szAtomName1[12], szAtomName2[12];
01095 int nAtomSiteId1, nAtomSiteId2;
01096 int nComponentId1, nComponentId2;
01097 int nChainId1, nChainId2;
01098 int nActiveSiteId1, nActiveSiteId2, nActiveSiteId3;
01099 int nAssemblyId, nAssemblySerial;
01100 int nModelId;
01101 char szCode3letter1[4];
01102 char szCode3letter2[4];
01103 char szChainID1[4];
01104 char szChainID2[4];
01105 int nPdbSeq1, nPdbSeq2;
01106 int nCarbonIndex = getElementNo(" C");
01107 int nNeighbourTypeId;
01108 int nNeighbourTypeId1;
01109 int nNeighbourTypeId2;
01110 int nOrdering;
01111 char chLigandAcceptorDonor[2];
01112 char chNeighbourAcceptorDonor[2];
01113 char szSymbol1[3];
01114 char szSymbol2[3];
01115
01116
01117 cout << " ------------- start " << nDepId << " ---------------------" << endl;
01118 MSDManager mgr;
01119 mgr.LoadWarehouse(db, nDepId);
01120 cout << " ------------- loaded " << nDepId << " ---------------------" << endl;
01121
01122
01123
01124 sprintf(select_stm,"select /*+ INDEX(ENTRY)*/ entry.SHORT_EXPERIMENT_TYPE "
01125 " from entry where entry.entry_id = :nDepDepId <int> ");
01126 double dblMinDist;
01127 int nSeqDist;
01128 try{
01129 msdbStream imsd (50, select_stm, db);
01130 imsd << nDepDepId;
01131
01132 while (!imsd.eof()){
01133 imsd >> szExpType;
01134 nExpTypeInd = imsd.is_null();
01135 if (!nExpTypeInd)
01136 trimRight(szExpType);
01137 else
01138 szExpType[0] = 0;
01139 ExpType::EExpType nType = ExpType::nUNKNOWN;
01140 if (!strcmp("XRAY", szExpType))
01141 nType = ExpType::nXRAY;
01142 else if (!strcmp("NMR", szExpType))
01143 nType = ExpType::nNMR;
01144 else if (!strcmp("SSNMR", szExpType))
01145 nType = ExpType::nNMR;
01146 int nCells = nType == ExpType::nXRAY? nMaxCells : -1;
01147 if (nCells >= 0 && mgr.IsTheoretical() == S_OK)
01148 nCells = -1;
01149
01150 dblMinDist = 0.1;
01151 nSeqDist = 1;
01152
01153 if (nCells >= 0)
01154 {
01155 if (!mgr.isCrystInfo())
01156 return Error("Missing cell parameters (a,b,c,alpha,beta,gamma)", E_FAIL);
01157 else if (!mgr.isSpaceGroup())
01158 return Error("Missing space group of symmetry", E_FAIL);
01159 else if (!mgr.isTransfMatrix())
01160 return Error("Orthogonalization/fractionalization is impossible", E_FAIL);
01161 }
01162 }
01163 }
01164 catch(msdbException& p){
01165 cerr<<"get experiment type error"<<endl;
01166 cerr<<p.msg<<endl;
01167 cerr<<p.stm_text<<endl;
01168 cerr<<p.var_info<<endl;
01169 }
01170
01171 CDBAssembly** vpAssembly = NULL;
01172 mgr.GetAssemblyArray(&vpAssembly);
01173 int nAssemblyCount = 0;
01174 mgr.GetAssemblyCount(&nAssemblyCount);
01175 CDBAtom** vpMacroMolAtom = new CDBAtom*[mgr.GetAtomArrayLength()];
01176 CDBAtom** vpLigandAtom = new CDBAtom*[mgr.GetAtomArrayLength()];
01177 CDBAtom** vpWaterAtom = new CDBAtom*[mgr.GetAtomArrayLength()];
01178 CAtom** vpMaxOccAtom = new CAtom*[mgr.GetAtomArrayLength()];
01179 CAtom** vpFirstReserveAtom = new CAtom*[mgr.GetAtomArrayLength()];
01180 CAtom** vpSecondReserveAtom = new CAtom*[mgr.GetAtomArrayLength()];
01181 int* vnActiveSiteId = new int[mgr.GetAtomArrayLength() + 1];
01182 for (int l = 0; l < nAssemblyCount; l++)
01183 {
01184 CDBAssembly* pAssembly = vpAssembly[l];
01185 pAssembly->GetId(&nAssemblyId);
01186 pAssembly->GetSerial(&nAssemblySerial);
01187 int nModelAssemblyCount = 0;
01188 pAssembly->GetModelAssemblyCount(&nModelAssemblyCount);
01189 for (int im = 0; im < nModelAssemblyCount; im++)
01190 {
01191 memset(vnActiveSiteId, 0, sizeof(int) * (mgr.GetAtomArrayLength() + 1));
01192 CDBModelAssembly* pModelAssembly = NULL;
01193 pAssembly->GetModelAssembly(&pModelAssembly, im);
01194 CDBModel* pModel = NULL;
01195 pModelAssembly->GetModel(&pModel);
01196 pModel->GetId(&nModelId);
01197 int nMacroMolAtoms = 0;
01198 int nLigandAtoms = 0;
01199 int nWaterAtoms = 0;
01200 int nMaxOccAtoms = 0;
01201 for (int nc = 0; nc < pModel->GetNumberOfChains(); nc++)
01202 {
01203 CDBChain* pChain = (CDBChain*)pModel->GetChain(nc);
01204 char szType[20];
01205 pChain->GetType(szType);
01206 bool bIsWater = (szType[0] != 0 && szType[1] == 'W') || szType[0] == 'W';
01207 bool bIsLigand = (szType[0] != 0 && szType[1] == 'B') || szType[0] == 'B';
01208 bool bIsMacromol = (szType[0] != 0 && szType[1] == 'C') || szType[0] == 'C';
01209 bool bIsModRes = false;
01210 if (bIsLigand || bIsMacromol || bIsWater) for (int nr = 0; nr < pChain->GetNumberOfResidues(); nr++)
01211 {
01212 CDBResidue* pResidue = (CDBResidue*)pChain->GetResidue(nr);
01213 if (bIsMacromol && pChain->IsProtein() == S_OK)
01214 bIsModRes = !isStdAminoAcid(pResidue->name);
01215 int nResAssemblyId = 0;
01216 if (pResidue->GetAssemblyId(&nResAssemblyId) != S_OK ||
01217 nAssemblyId != nResAssemblyId)
01218 continue;
01219 rvector occupancy = NULL;
01220 int nAltLocs,alflag;
01221 AltLoc aLoc;
01222 PAltLoc aL = NULL;
01223 pResidue->GetAltLocations ( nAltLocs, aL, occupancy, alflag );
01224 PCAtom pPrevAtom = NULL;
01225 double dx = -2.0;
01226 aLoc[0] = char(0);
01227 if (nAltLocs > 1)
01228 {
01229 for (int i=0;i<nAltLocs;i++)
01230 if ((aL[i][0]) && (occupancy[i]>dx)) {
01231 dx = occupancy[i];
01232 strcpy ( aLoc,aL[i] );
01233 }
01234 }
01235 else
01236 {
01237 if (strcmp(aLoc,aL[0]))
01238 strcpy ( aLoc,aL[0] );
01239 }
01240 for (int na = 0; na < pResidue->nAtoms; na++)
01241 {
01242 CDBAtom* pAtom = (CDBAtom*)pResidue->atom[na];
01243 if (!pAtom || pAtom->Ter)
01244 continue;
01245 if (bIsLigand || bIsModRes)
01246 vpLigandAtom[nLigandAtoms++] = pAtom;
01247 else if (bIsMacromol && !bIsModRes)
01248 vpMacroMolAtom[nMacroMolAtoms++] = pAtom;
01249 else if (bIsWater)
01250 vpWaterAtom[nWaterAtoms++] = pAtom;
01251 if (pAtom->element[0] == 0)
01252 continue;
01253 if (nAltLocs > 1)
01254 {
01255 bool B = !strcmp(aLoc, pAtom->altLoc);
01256 if ((!B) && (!pAtom->altLoc[0])) {
01257
01258
01259 for (int ja=na+1;(ja<pResidue->nAtoms) && (!B);ja++)
01260 if (pResidue->atom[ja]) {
01261 if ((!pResidue->atom[ja]->Ter) &&
01262 (!strcmp(pResidue->atom[ja]->name, pAtom->name)))
01263 B = !strcmp(aLoc,pResidue->atom[ja]->altLoc);
01264 }
01265
01266
01267 B = !B;
01268 }
01269 if (B)
01270 vpMaxOccAtom[nMaxOccAtoms++] = pAtom;
01271 }
01272 else
01273 vpMaxOccAtom[nMaxOccAtoms++] = pAtom;
01274 }
01275 }
01276 }
01277
01278 if (!nMacroMolAtoms || !nLigandAtoms)
01279 continue;
01280
01281
01282 donor_out* d_out = new donor_out[nMaxOccAtoms];
01283 next_d* nextd = new next_d[nMaxOccAtoms];
01284 int nDonorCount = 0;
01285 acceptor_out* a_out = new acceptor_out[nMaxOccAtoms];
01286 next_a* nexta = new next_a[nMaxOccAtoms];
01287 int nAcceptorCount = 0;
01288 std::map<CAtomKey*, int, ltatom> mapDonor;
01289 initDonorOut(vpMaxOccAtom, nMaxOccAtoms,
01290 connect_d, nDCount, mapDonor,
01291 d_out, nextd, &nDonorCount);
01292 std::map<CAtomKey*, int, ltatom> mapAcceptor;
01293 initAcceptorOut(vpMaxOccAtom, nMaxOccAtoms,
01294 connect_a, nACount, mapAcceptor,
01295 a_out, nexta, &nAcceptorCount);
01296 PSContact vContact = NULL;
01297 int nContacts = 0;
01298
01299 mgr.CMMDBManager::SeekContacts((PPCAtom)vpLigandAtom, nLigandAtoms,
01300 (PPCAtom)vpMacroMolAtom, nMacroMolAtoms,
01301 dblMinDist, dblMaxDist, nSeqDist,
01302 vContact, nContacts, 0, NULL, 1);
01303 mgr.CMMDBManager::SeekContacts((PPCAtom)vpLigandAtom, nLigandAtoms,
01304 (PPCAtom)vpLigandAtom, nLigandAtoms,
01305 dblMinDist, dblMaxDist, nSeqDist,
01306 vContact, nContacts, 0, NULL, 2);
01307 mgr.CMMDBManager::SeekContacts((PPCAtom)vpLigandAtom, nLigandAtoms,
01308 (PPCAtom)vpWaterAtom, nWaterAtoms,
01309 dblMinDist, dblMaxDist, nSeqDist,
01310 vContact, nContacts, 0, NULL, 3);
01311 int *vnContactType = nContacts > 0 ? new int[nContacts] : NULL;
01312 if (vnContactType != NULL)
01313 memset(vnContactType, 0, sizeof(int) * nContacts);
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337 std::map<int, double> ligandDist;
01338 for (int ic = 0, nCount = 0; ic < nContacts; ic ++) if (vContact[ic].group != 3)
01339 {
01340 int nLigandAtom = vContact[ic].id1;
01341 int nNeighbourAtom = vContact[ic].id2;
01342 CDBAtom* pLigandAtom = vpLigandAtom[nLigandAtom];
01343 CDBResidue* pLigandResidue = (CDBResidue*)pLigandAtom->residue;
01344 int nId = 0;
01345 pLigandResidue->GetId(&nId);
01346 std::map<int, double>::iterator it = ligandDist.find(nId);
01347 if (it != ligandDist.end() && (*it).second < 4.)
01348 continue;
01349 ligandDist[nId] = vContact[ic].dist;
01350 }
01351 for (int ic = 0, nCount = 0; ic < nContacts; ic ++)
01352 {
01353
01354 int nLigandAtom = vContact[ic].id1;
01355 CDBAtom* pLigandAtom = vpLigandAtom[nLigandAtom];
01356 CDBResidue* pLigandResidue = (CDBResidue*)pLigandAtom->residue;
01357 CDBChain* pLigandChain = (CDBChain*)pLigandResidue->chain;
01358 int nLigandResidueId = 0;
01359 pLigandResidue->GetId(&nLigandResidueId);
01360
01361 int nNeighbourAtom = vContact[ic].id2;
01362 CDBAtom* pNeighbourAtom = vContact[ic].group == 1
01363 ? vpMacroMolAtom[nNeighbourAtom]
01364 : vContact[ic].group == 2 ? vpLigandAtom[nNeighbourAtom]
01365 : vpWaterAtom[nNeighbourAtom];
01366 CDBResidue* pNeighbourResidue = (CDBResidue*)pNeighbourAtom->residue;
01367
01368 bool bExtendedDistance = false;
01369 if (vContact[ic].group != 3)
01370 {
01371 std::map<int, double>::iterator itLigand = ligandDist.find(nLigandResidueId);
01372 bExtendedDistance = (*itLigand).second >= 4.;
01373 }
01374
01375 int nAtomNo1 = getElementNo(pLigandAtom->element);
01376 int nAtomNo2 = getElementNo(pNeighbourAtom->element);
01377 EContactType nBondType = UNDEFINED_BOND;
01378 bool bLigandIsDonor = false;
01379 if (!(nAtomNo1 == ELEMENT_UNKNOWN || nAtomNo2 == ELEMENT_UNKNOWN))
01380 {
01381 nBondType = getBondType(pLigandAtom, pNeighbourAtom, vContact[ic].dist);
01382 if (nBondType == PROHIBITED_BOND)
01383 continue;
01384 if ((pLigandAtom->residue->nAtoms > 1 || !strcmp(pLigandAtom->residue->name, "HOH")) &&
01385 (pNeighbourAtom->residue->nAtoms > 1 || !strcmp(pNeighbourAtom->residue->name, "HOH")) &&
01386 (nBondType == UNDEFINED_BOND || nBondType == VAN_DEF_WAALS_BOND) && vContact[ic].dist < 4.)
01387 {
01388 bool bHydrogen = false;
01389 bool bHydrogen2= false;
01390 bool bDisulfid = false;
01391 CAtom** vpNeighbourAtom = NULL;
01392 int nNeighbourAtoms = 0;
01393 switch(vContact[ic].group)
01394 {
01395 case 1:
01396 vpNeighbourAtom = (CAtom**)vpMacroMolAtom;
01397 nNeighbourAtoms = nMacroMolAtoms;
01398 break;
01399 case 2:
01400 {
01401 int nSourceArraySize = nLigandAtoms;
01402 CAtom** vpSourceArray = (CAtom**)vpLigandAtom;
01403 vpNeighbourAtom = vpFirstReserveAtom;
01404 int nStartAtom = nNeighbourAtom - 1;
01405 int nNeighbourIndex = ((CAtomWrapper*)pNeighbourAtom)->GetIndex();
01406 while (nStartAtom >= 0 &&
01407 vpSourceArray[nStartAtom]->residue->chain == pNeighbourAtom->residue->chain &&
01408 ((CAtomWrapper*)vpSourceArray[nStartAtom])->GetIndex() < nNeighbourIndex)
01409 nStartAtom--;
01410 int nEndAtom = nNeighbourAtom + 1;
01411 while (nEndAtom < nSourceArraySize &&
01412 vpSourceArray[nEndAtom]->residue->chain == pNeighbourAtom->residue->chain &&
01413 ((CAtomWrapper*)vpSourceArray[nEndAtom])->GetIndex() > nNeighbourIndex)
01414 nEndAtom++;
01415 for (int iatm = nStartAtom + 1, nNeighbourAtoms = 0; iatm < nEndAtom; iatm++)
01416 vpNeighbourAtom[nNeighbourAtoms++] = vpSourceArray[iatm];
01417 break;
01418 }
01419 case 3:
01420 {
01421 int nSourceArraySize = nWaterAtoms;
01422 CAtom** vpSourceArray = (CAtom**)vpWaterAtom;
01423 vpNeighbourAtom = vpFirstReserveAtom;
01424 int nStartAtom = nNeighbourAtom - 1;
01425 int nNeighbourIndex = ((CAtomWrapper*)pNeighbourAtom)->GetIndex();
01426 while (nStartAtom >= 0 &&
01427 vpSourceArray[nStartAtom]->residue == pNeighbourAtom->residue &&
01428 ((CAtomWrapper*)vpSourceArray[nStartAtom])->GetIndex() < nNeighbourIndex)
01429 nStartAtom--;
01430 int nEndAtom = nNeighbourAtom + 1;
01431 while (nEndAtom < nSourceArraySize &&
01432 vpSourceArray[nEndAtom]->residue == pNeighbourAtom->residue &&
01433 ((CAtomWrapper*)vpSourceArray[nEndAtom])->GetIndex() > nNeighbourIndex)
01434 nEndAtom++;
01435 for (int iatm = nStartAtom + 1, nNeighbourAtoms = 0; iatm < nEndAtom; iatm++)
01436 vpNeighbourAtom[nNeighbourAtoms++] = vpSourceArray[iatm];
01437 }
01438 break;
01439 }
01440
01441 int nLigandChainAtomCount = 0;
01442 PPCAtom vpLigandChainAtom = vpSecondReserveAtom;
01443 int nStartAtom = nLigandAtom - 1;
01444 int nLigandIndex = ((CAtomWrapper*)pLigandAtom)->GetIndex();
01445 while (nStartAtom >= 0 &&
01446 vpLigandAtom[nStartAtom]->residue->chain == pLigandResidue->chain &&
01447 ((CAtomWrapper*)vpLigandAtom[nStartAtom])->GetIndex() < nLigandIndex)
01448 nStartAtom--;
01449 int nEndAtom = nLigandAtom + 1;
01450 while (nEndAtom < nLigandAtoms &&
01451 vpLigandAtom[nEndAtom]->residue->chain == pLigandResidue->chain &&
01452 ((CAtomWrapper*)vpLigandAtom[nEndAtom])->GetIndex() > nLigandIndex)
01453 nEndAtom++;
01454 for (int iatm = nStartAtom + 1, nLigandChainAtomCount = 0; iatm < nEndAtom; iatm++)
01455 vpLigandChainAtom[nLigandChainAtomCount++] = vpLigandAtom[iatm];
01456 int nLigandDonorNumber = -1;
01457 donor_out ligandDonor;
01458 next_d ligandNextD;
01459 Atom_GetDonorInfo(pLigandAtom, (CAtom**)vpMaxOccAtom, mapDonor, d_out, nextd, nDonorCount,
01460 vpLigandChainAtom, nLigandChainAtomCount,
01461 nLigandDonorNumber, ligandDonor, ligandNextD);
01462 int nLigandAcceptorNumber = -1;
01463 acceptor_out ligandAcceptor;
01464 next_a ligandNextA;
01465 Atom_GetAcceptorInfo(pLigandAtom, (CAtom**)vpMaxOccAtom, mapAcceptor, a_out, nexta, nAcceptorCount,
01466 vpLigandChainAtom, nLigandChainAtomCount,
01467 nLigandAcceptorNumber, ligandAcceptor, ligandNextA);
01468 int nMacroDonorNumber = -1;
01469 donor_out macroDonor;
01470 next_d macroNextD;
01471 Atom_GetDonorInfo(pNeighbourAtom, (CAtom**)vpMaxOccAtom, mapDonor, d_out, nextd, nDonorCount,
01472 vpNeighbourAtom, nNeighbourAtoms,
01473 nMacroDonorNumber, macroDonor, macroNextD);
01474 int nMacroAcceptorNumber = -1;
01475 acceptor_out macroAcceptor;
01476 next_a macroNextA;
01477 Atom_GetAcceptorInfo(pNeighbourAtom, (CAtom**)vpMaxOccAtom, mapAcceptor, a_out, nexta, nAcceptorCount,
01478 vpNeighbourAtom, nNeighbourAtoms,
01479 nMacroAcceptorNumber, macroAcceptor, macroNextA);
01480 if ((nMacroDonorNumber == -1 || nLigandAcceptorNumber == -1) &&
01481 (nMacroAcceptorNumber == -1 || nLigandDonorNumber == -1))
01482 {
01483
01484 if (nLigandAcceptorNumber == -1 && nLigandDonorNumber == -1)
01485 {
01486 bool bAtomXIsDonor = false;
01487 bool bAtomXIsAcceptor = false;
01488 PCCalcGraph pGraph = NULL;
01489 int nSupChemCompId = 0;
01490 pLigandResidue->GetSupChemCompId(&nSupChemCompId);
01491 mgr.GetSupChemCompGraph(nSupChemCompId, &pGraph);
01492 if (pGraph)
01493 Atom_GetDonorAcceptorInfo(pLigandAtom,
01494 2,
01495 pGraph,
01496 nAtomNo1,
01497 (CDBAtom**)vpLigandChainAtom,
01498 nLigandChainAtomCount,
01499 ligandDonor,
01500 ligandNextD,
01501 ligandAcceptor,
01502 ligandNextA,
01503 bAtomXIsDonor,
01504 bAtomXIsAcceptor);
01505 if (bAtomXIsDonor)
01506 nLigandDonorNumber = 0;
01507 if (bAtomXIsAcceptor)
01508 nLigandAcceptorNumber = 0;
01509 }
01510 if (nMacroAcceptorNumber == -1 && nMacroDonorNumber == -1)
01511 {
01512 bool bAtomXIsDonor = false;
01513 bool bAtomXIsAcceptor = false;
01514 PCCalcGraph pGraph = NULL;
01515 int nSupChemCompId = 0;
01516 pNeighbourResidue->GetSupChemCompId(&nSupChemCompId);
01517 mgr.GetSupChemCompGraph(nSupChemCompId, &pGraph);
01518 if (pGraph)
01519 Atom_GetDonorAcceptorInfo(pNeighbourAtom,
01520 vContact[ic].group,
01521 pGraph,
01522 nAtomNo2,
01523 (CDBAtom**)vpNeighbourAtom,
01524 nNeighbourAtoms,
01525 macroDonor,
01526 macroNextD,
01527 macroAcceptor,
01528 macroNextA,
01529 bAtomXIsDonor,
01530 bAtomXIsAcceptor);
01531 if (bAtomXIsDonor)
01532 nMacroDonorNumber = 0;
01533 if (bAtomXIsAcceptor)
01534 nMacroAcceptorNumber = 0;
01535 }
01536 }
01537 if (nMacroDonorNumber != -1 && nLigandAcceptorNumber != -1)
01538 {
01539 bHydrogen = isHBond(macroDonor, macroNextD,
01540 ligandAcceptor, ligandNextA,
01541 vContact[ic].dist, bDisulfid);
01542 bLigandIsDonor = false;
01543 }
01544 if (!bHydrogen && !bDisulfid && nMacroAcceptorNumber != -1 && nLigandDonorNumber != -1)
01545 {
01546 bHydrogen = isHBond(ligandDonor, ligandNextD,
01547 macroAcceptor, macroNextA,
01548 vContact[ic].dist, bDisulfid);
01549 bLigandIsDonor = true;
01550 }
01551 if (nMacroDonorNumber == -1 && nLigandAcceptorNumber == -1 ||
01552 nMacroAcceptorNumber == -1 && nLigandDonorNumber == -1)
01553 {
01554 bHydrogen2 = isHydrogenBond(pLigandAtom, pNeighbourAtom, vContact[ic].dist, &bLigandIsDonor);
01555 }
01556 if (bHydrogen)
01557 nBondType = HYDROGEN_BOND_HBEXPLORE;
01558 else if (bHydrogen2)
01559 nBondType = HYDROGEN_BOND_BNDLIST;
01560 else if (bDisulfid)
01561 nBondType = DISULFIDE_BOND;
01562 }
01563 }
01564 nActiveSiteType1 =
01565 nAtomNo1 == ELEMENT_UNKNOWN &&
01566 !strcmp("GVC", pLigandAtom->GetAtomName())
01567 ? 2 : 1;
01568 nActiveSiteType2 =
01569 nAtomNo2 == ELEMENT_UNKNOWN &&
01570 !strcmp("GVC", pNeighbourAtom->GetAtomName())
01571 ? 2 : 1;
01572 if (nActiveSiteType1 != nActiveSiteType2)
01573 continue;
01574 if(nActiveSiteType1 == 2)
01575 nBondType = TWO_PLANES_BOND;
01576 if (nBondType != TWO_PLANES_BOND && vContact[ic].dist > 4. && !bExtendedDistance)
01577 continue;
01578
01579
01580
01581
01582
01583 }
01584
01585
01586
01587
01588
01589
01590 if (vnContactType)
01591 delete vnContactType;
01592 if (vContact)
01593 delete vContact;
01594 clear(mapDonor);
01595 clear(mapAcceptor);
01596 }
01597 }
01598
01599 delete [] vnActiveSiteId;
01600 delete [] vpWaterAtom;
01601 delete [] vpMacroMolAtom;
01602 delete [] vpLigandAtom;
01603 delete [] vpFirstReserveAtom;
01604 delete [] vpSecondReserveAtom;
01605 return S_OK;
01606 }
01607
01608 #endif
01609
01610
