atom_set.cpp
上传用户:yhdzpy8989
上传日期:2007-06-13
资源大小:13604k
文件大小:15k
源码类别:

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: atom_set.cpp,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 18:27:47  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.20
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: atom_set.cpp,v 1000.2 2004/06/01 18:27:47 gouriano Exp $
  10. * ===========================================================================
  11. *
  12. *                            PUBLIC DOMAIN NOTICE
  13. *               National Center for Biotechnology Information
  14. *
  15. *  This software/database is a "United States Government Work" under the
  16. *  terms of the United States Copyright Act.  It was written as part of
  17. *  the author's official duties as a United States Government employee and
  18. *  thus cannot be copyrighted.  This software/database is freely available
  19. *  to the public for use. The National Library of Medicine and the U.S.
  20. *  Government have not placed any restriction on its use or reproduction.
  21. *
  22. *  Although all reasonable efforts have been taken to ensure the accuracy
  23. *  and reliability of the software and data, the NLM and the U.S.
  24. *  Government do not and cannot warrant the performance or results that
  25. *  may be obtained by using this software or data. The NLM and the U.S.
  26. *  Government disclaim all warranties, express or implied, including
  27. *  warranties of performance, merchantability or fitness for any particular
  28. *  purpose.
  29. *
  30. *  Please cite the author in any work or product based on this material.
  31. *
  32. * ===========================================================================
  33. *
  34. * Authors:  Paul Thiessen
  35. *
  36. * File Description:
  37. *      Classes to hold sets of atomic data
  38. *
  39. * ===========================================================================
  40. */
  41. #ifdef _MSC_VER
  42. #pragma warning(disable:4018)   // disable signed/unsigned mismatch warning in MSVC
  43. #endif
  44. #include <ncbi_pch.hpp>
  45. #include <corelib/ncbistd.hpp>
  46. #include <objects/mmdb1/Molecule_id.hpp>
  47. #include <objects/mmdb1/Residue_id.hpp>
  48. #include <objects/mmdb1/Atom_id.hpp>
  49. #include <objects/mmdb2/Model_space_points.hpp>
  50. #include <objects/mmdb2/Atomic_temperature_factors.hpp>
  51. #include <objects/mmdb2/Atomic_occupancies.hpp>
  52. #include <objects/mmdb2/Alternate_conformation_ids.hpp>
  53. #include <objects/mmdb2/Alternate_conformation_id.hpp>
  54. #include <objects/mmdb2/Anisotro_temperatu_factors.hpp>
  55. #include <objects/mmdb2/Isotropi_temperatu_factors.hpp>
  56. #include <objects/mmdb2/Conformation_ensemble.hpp>
  57. #include <objects/mmdb3/Atom_pntrs.hpp>
  58. #include "atom_set.hpp"
  59. #include "vector_math.hpp"
  60. #include "cn3d_tools.hpp"
  61. #include "structure_set.hpp"
  62. USING_NCBI_SCOPE;
  63. USING_SCOPE(objects);
  64. BEGIN_SCOPE(Cn3D)
  65. AtomSet::AtomSet(StructureBase *parent, const CAtomic_coordinates& coords) :
  66.     StructureBase(parent), activeEnsemble(NULL)
  67. {
  68.     int nAtoms = coords.GetNumber_of_points();
  69.     TRACEMSG("model has " << nAtoms << " atomic coords");
  70.     bool haveTemp = coords.IsSetTemperature_factors(),
  71.          haveOccup = coords.IsSetOccupancies(),
  72.          haveAlt = coords.IsSetAlternate_conf_ids();
  73.     // sanity check
  74.     if (coords.GetAtoms().GetMolecule_ids().size()!=nAtoms ||
  75.         coords.GetAtoms().GetResidue_ids().size()!=nAtoms ||
  76.         coords.GetAtoms().GetAtom_ids().size()!=nAtoms ||
  77.         coords.GetSites().GetX().size()!=nAtoms ||
  78.         coords.GetSites().GetY().size()!=nAtoms ||
  79.         coords.GetSites().GetZ().size()!=nAtoms ||
  80.         (haveTemp &&
  81.             ((coords.GetTemperature_factors().IsIsotropic() &&
  82.               coords.GetTemperature_factors().GetIsotropic().GetB().size()!=nAtoms) ||
  83.              (coords.GetTemperature_factors().IsAnisotropic() &&
  84.               (coords.GetTemperature_factors().GetAnisotropic().GetB_11().size()!=nAtoms ||
  85.                coords.GetTemperature_factors().GetAnisotropic().GetB_12().size()!=nAtoms ||
  86.                coords.GetTemperature_factors().GetAnisotropic().GetB_13().size()!=nAtoms ||
  87.                coords.GetTemperature_factors().GetAnisotropic().GetB_22().size()!=nAtoms ||
  88.                coords.GetTemperature_factors().GetAnisotropic().GetB_23().size()!=nAtoms ||
  89.                coords.GetTemperature_factors().GetAnisotropic().GetB_33().size()!=nAtoms)))) ||
  90.         (haveOccup && coords.GetOccupancies().GetO().size()!=nAtoms) ||
  91.         (haveAlt && coords.GetAlternate_conf_ids().Get().size()!=nAtoms))
  92.         ERRORMSG("AtomSet: confused by list length mismatch");
  93.     // atom-pntr iterators
  94.     CAtom_pntrs::TMolecule_ids::const_iterator
  95.         i_mID = coords.GetAtoms().GetMolecule_ids().begin();
  96.     CAtom_pntrs::TResidue_ids::const_iterator
  97.         i_rID = coords.GetAtoms().GetResidue_ids().begin();
  98.     CAtom_pntrs::TAtom_ids::const_iterator
  99.         i_aID = coords.GetAtoms().GetAtom_ids().begin();
  100.     // atom site iterators
  101.     double siteScale = static_cast<double>(coords.GetSites().GetScale_factor());
  102.     CModel_space_points::TX::const_iterator i_X = coords.GetSites().GetX().begin();
  103.     CModel_space_points::TY::const_iterator i_Y = coords.GetSites().GetY().begin();
  104.     CModel_space_points::TZ::const_iterator i_Z = coords.GetSites().GetZ().begin();
  105.     // occupancy iterator
  106.     CAtomic_occupancies::TO::const_iterator i_occup;
  107.     double occupScale = 0.0;
  108.     if (haveOccup) {
  109.         occupScale = static_cast<double>(coords.GetOccupancies().GetScale_factor());
  110.         i_occup = coords.GetOccupancies().GetO().begin();
  111.     }
  112.     // altConfID iterator
  113.     CAlternate_conformation_ids::Tdata::const_iterator i_alt;
  114.     if (haveAlt) i_alt = coords.GetAlternate_conf_ids().Get().begin();
  115.     // temperature iterators
  116.     double tempScale = 0.0;
  117.     CIsotropic_temperature_factors::TB::const_iterator i_tempI;
  118.     CAnisotropic_temperature_factors::TB_11::const_iterator i_tempA11;
  119.     CAnisotropic_temperature_factors::TB_12::const_iterator i_tempA12;
  120.     CAnisotropic_temperature_factors::TB_13::const_iterator i_tempA13;
  121.     CAnisotropic_temperature_factors::TB_22::const_iterator i_tempA22;
  122.     CAnisotropic_temperature_factors::TB_23::const_iterator i_tempA23;
  123.     CAnisotropic_temperature_factors::TB_33::const_iterator i_tempA33;
  124.     if (haveTemp) {
  125.         if (coords.GetTemperature_factors().IsIsotropic()) {
  126.             tempScale = static_cast<double>
  127.                 (coords.GetTemperature_factors().GetIsotropic().GetScale_factor());
  128.             i_tempI = coords.GetTemperature_factors().GetIsotropic().GetB().begin();
  129.         } else {
  130.             tempScale = static_cast<double>
  131.                 (coords.GetTemperature_factors().GetAnisotropic().GetScale_factor());
  132.             i_tempA11 = coords.GetTemperature_factors().GetAnisotropic().GetB_11().begin();
  133.             i_tempA12 = coords.GetTemperature_factors().GetAnisotropic().GetB_12().begin();
  134.             i_tempA13 = coords.GetTemperature_factors().GetAnisotropic().GetB_13().begin();
  135.             i_tempA22 = coords.GetTemperature_factors().GetAnisotropic().GetB_22().begin();
  136.             i_tempA23 = coords.GetTemperature_factors().GetAnisotropic().GetB_23().begin();
  137.             i_tempA33 = coords.GetTemperature_factors().GetAnisotropic().GetB_33().begin();
  138.         }
  139.     }
  140.     const StructureObject *constObject;
  141.     if (!GetParentOfType(&constObject)) return;
  142.     StructureObject *object = const_cast<StructureObject*>(constObject);
  143.     // actually do the work of unpacking serial atom data into Atom objects
  144.     for (int i=0; i<nAtoms; ++i) {
  145.         AtomCoord *atom = new AtomCoord(this);
  146.         atom->site.x = (static_cast<double>(*(i_X++)))/siteScale;
  147.         atom->site.y = (static_cast<double>(*(i_Y++)))/siteScale;
  148.         atom->site.z = (static_cast<double>(*(i_Z++)))/siteScale;
  149.         if (haveOccup)
  150.             atom->occupancy = (static_cast<double>(*(i_occup++)))/occupScale;
  151.         if (haveAlt)
  152.             atom->altConfID = (i_alt++)->Get()[0];
  153.         if (haveTemp) {
  154.             if (coords.GetTemperature_factors().IsIsotropic()) {
  155.                 atom->averageTemperature =
  156.                     (static_cast<double>(*(i_tempI++)))/tempScale;
  157.             } else {
  158.                 atom->averageTemperature = (
  159.                     (static_cast<double>(*(i_tempA11++))) +
  160.                     (static_cast<double>(*(i_tempA12++))) +
  161.                     (static_cast<double>(*(i_tempA13++))) +
  162.                     (static_cast<double>(*(i_tempA22++))) +
  163.                     (static_cast<double>(*(i_tempA23++))) +
  164.                     (static_cast<double>(*(i_tempA33++)))) / (tempScale * 6.0);
  165.             }
  166.             // track min and max temperatures over whole object
  167.             if (object->minTemperature == StructureObject::NO_TEMPERATURE ||
  168.                 atom->averageTemperature < object->minTemperature)
  169.                 object->minTemperature = atom->averageTemperature;
  170.             if (object->maxTemperature == StructureObject::NO_TEMPERATURE ||
  171.                 atom->averageTemperature > object->maxTemperature)
  172.                 object->maxTemperature = atom->averageTemperature;
  173.         }
  174.         // store pointer in map - key+altConfID must be unique
  175.         const AtomPntrKey& key = MakeKey(AtomPntr(
  176.             *(i_mID++),
  177.             *(i_rID++),
  178.             *(i_aID++)));
  179.         if (atomMap.find(key) != atomMap.end()) {
  180.             AtomAltList::const_iterator i_atom, e=atomMap[key].end();
  181.             for (i_atom=atomMap[key].begin(); i_atom!=e; ++i_atom) {
  182.                 if ((*i_atom)->altConfID == atom->altConfID)
  183.                     ERRORMSG("confused by multiple atoms of same pntr+altConfID");
  184.             }
  185.         }
  186.         atomMap[key].push_back(atom);
  187.         if (i==0) TRACEMSG("first atom: x " << atom->site.x <<
  188.                 ", y " << atom->site.y <<
  189.                 ", z " << atom->site.z <<
  190.                 ", occup " << atom->occupancy <<
  191.                 ", altConfId '" << atom->altConfID << "'" <<
  192.                 ", temp " << atom->averageTemperature);
  193.     }
  194.     // get alternate conformer ensembles; store as string of altID characters
  195.     if (haveAlt && coords.IsSetConf_ensembles()) {
  196.         CAtomic_coordinates::TConf_ensembles::const_iterator i_ensemble,
  197.             e_ensemble = coords.GetConf_ensembles().end();
  198.         for (i_ensemble=coords.GetConf_ensembles().begin(); i_ensemble!=e_ensemble; ++i_ensemble) {
  199.             const CConformation_ensemble& ensemble = (*i_ensemble).GetObject();
  200.             string *ensembleStr = new string();
  201.             CConformation_ensemble::TAlt_conf_ids::const_iterator i_altIDs,
  202.                 e_altIDs = ensemble.GetAlt_conf_ids().end();
  203.             for (i_altIDs=ensemble.GetAlt_conf_ids().begin();
  204.                  i_altIDs!=e_altIDs; ++i_altIDs) {
  205.                 (*ensembleStr) += i_altIDs->Get()[0];
  206.             }
  207.             ensembles.push_back(ensembleStr);
  208.             TRACEMSG("alt conf ensemble '" << (*ensembleStr) << "'");
  209.         }
  210.     }
  211. }
  212. AtomSet::~AtomSet(void)
  213. {
  214.     EnsembleList::iterator i, e=ensembles.end();
  215.     for (i=ensembles.begin(); i!=e; ++i)
  216.         delete const_cast<string *>(*i);
  217. }
  218. const double AtomCoord::NO_TEMPERATURE = -1.0;
  219. const double AtomCoord::NO_OCCUPANCY = -1.0;
  220. const char AtomCoord::NO_ALTCONFID = '-';
  221. bool AtomSet::SetActiveEnsemble(const string *ensemble)
  222. {
  223.     // if not NULL, make sure it's one of this AtomSet's ensembles
  224.     if (ensemble) {
  225.         EnsembleList::const_iterator e, ee=ensembles.end();
  226.         for (e=ensembles.begin(); e!=ee; ++e) {
  227.             if (*e == ensemble) break;
  228.         }
  229.         if (e == ee) {
  230.             ERRORMSG("AtomSet::SetActiveEnsemble received invalid ensemble");
  231.             return false;
  232.         }
  233.     }
  234.     activeEnsemble = ensemble;
  235.     return true;
  236. }
  237. const AtomCoord* AtomSet::GetAtom(const AtomPntr& ap,
  238.     bool getAny, bool suppressWarning) const
  239. {
  240.     AtomMap::const_iterator atomConfs = atomMap.find(MakeKey(ap));
  241.     if (atomConfs == atomMap.end()) {
  242.         if (!suppressWarning)
  243.             WARNINGMSG("can't find atom(s) from pointer (" << ap.mID << ','
  244.                              << ap.rID << ',' << ap.aID << ')');
  245.         return NULL;
  246.     }
  247.     AtomAltList::const_iterator atom = atomConfs->second.begin();
  248.     // if no activeEnsemble specified, just return first altConf
  249.     if (!activeEnsemble) return *atom;
  250.     // otherwise, try to find first atom whose altConfID is in the activeEnsemble
  251.     AtomAltList::const_iterator e = atomConfs->second.end();
  252.     for (; atom!=e; ++atom) {
  253.         if (activeEnsemble->find((*atom)->altConfID) != activeEnsemble->npos)
  254.             return *atom;
  255.     }
  256.     // if none found, but getAny is true, just return the first of the list
  257.     if (getAny) return atomConfs->second.front();
  258.     return NULL;
  259. }
  260. AtomCoord::AtomCoord(StructureBase *parent) :
  261. StructureBase(parent),
  262.     averageTemperature(NO_TEMPERATURE),
  263.     occupancy(NO_OCCUPANCY),
  264.     altConfID(NO_ALTCONFID)
  265. {
  266. }
  267. END_SCOPE(Cn3D)
  268. /*
  269. * ---------------------------------------------------------------------------
  270. * $Log: atom_set.cpp,v $
  271. * Revision 1000.2  2004/06/01 18:27:47  gouriano
  272. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.20
  273. *
  274. * Revision 1.20  2004/05/21 21:41:38  gorelenk
  275. * Added PCH ncbi_pch.hpp
  276. *
  277. * Revision 1.19  2004/03/15 18:16:33  thiessen
  278. * prefer prefix vs. postfix ++/-- operators
  279. *
  280. * Revision 1.18  2004/02/19 17:04:42  thiessen
  281. * remove cn3d/ from include paths; add pragma to disable annoying msvc warning
  282. *
  283. * Revision 1.17  2003/10/21 13:48:48  grichenk
  284. * Redesigned type aliases in serialization library.
  285. * Fixed the code (removed CRef-s, added explicit
  286. * initializers etc.)
  287. *
  288. * Revision 1.16  2003/02/03 19:20:00  thiessen
  289. * format changes: move CVS Log to bottom of file, remove std:: from .cpp files, and use new diagnostic macros
  290. *
  291. * Revision 1.15  2001/12/12 14:04:12  thiessen
  292. * add missing object headers after object loader change
  293. *
  294. * Revision 1.14  2001/08/09 19:07:13  thiessen
  295. * add temperature and hydrophobicity coloring
  296. *
  297. * Revision 1.13  2001/06/02 17:22:45  thiessen
  298. * fixes for GCC
  299. *
  300. * Revision 1.12  2001/05/31 18:47:06  thiessen
  301. * add preliminary style dialog; remove LIST_TYPE; add thread single and delete all; misc tweaks
  302. *
  303. * Revision 1.11  2001/05/24 19:06:19  thiessen
  304. * fix to compile on SGI
  305. *
  306. * Revision 1.10  2000/08/25 14:22:00  thiessen
  307. * minor tweaks
  308. *
  309. * Revision 1.9  2000/08/03 15:12:22  thiessen
  310. * add skeleton of style and show/hide managers
  311. *
  312. * Revision 1.8  2000/07/27 13:30:51  thiessen
  313. * remove 'using namespace ...' from all headers
  314. *
  315. * Revision 1.7  2000/07/18 02:41:32  thiessen
  316. * fix bug in virtual bonds and altConfs
  317. *
  318. * Revision 1.6  2000/07/16 23:19:10  thiessen
  319. * redo of drawing system
  320. *
  321. * Revision 1.5  2000/07/12 23:27:48  thiessen
  322. * now draws basic CPK model
  323. *
  324. * Revision 1.4  2000/07/11 13:45:28  thiessen
  325. * add modules to parse chemical graph; many improvements
  326. *
  327. * Revision 1.3  2000/07/01 15:43:50  thiessen
  328. * major improvements to StructureBase functionality
  329. *
  330. * Revision 1.2  2000/06/29 19:17:47  thiessen
  331. * improved atom map
  332. *
  333. * Revision 1.1  2000/06/29 14:35:05  thiessen
  334. * new atom_set files
  335. *
  336. */