summaryrefslogtreecommitdiff
path: root/mmdb2/mmdb_coormngr.h
diff options
context:
space:
mode:
Diffstat (limited to 'mmdb2/mmdb_coormngr.h')
-rw-r--r--mmdb2/mmdb_coormngr.h959
1 files changed, 959 insertions, 0 deletions
diff --git a/mmdb2/mmdb_coormngr.h b/mmdb2/mmdb_coormngr.h
new file mode 100644
index 0000000..c493169
--- /dev/null
+++ b/mmdb2/mmdb_coormngr.h
@@ -0,0 +1,959 @@
+// $Id: mmdb_coormngr.h $
+// =================================================================
+//
+// CCP4 Coordinate Library: support of coordinate-related
+// functionality in protein crystallography applications.
+//
+// Copyright (C) Eugene Krissinel 2000-2013.
+//
+// This library is free software: you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License version 3, modified in accordance with the provisions
+// of the license to address the requirements of UK law.
+//
+// You should have received a copy of the modified GNU Lesser
+// General Public License along with this library. If not, copies
+// may be downloaded from http://www.ccp4.ac.uk/ccp4license.php
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU Lesser General Public License for more details.
+//
+// =================================================================
+//
+// 14.07.13 <-- Date of Last Modification.
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// -----------------------------------------------------------------
+//
+// **** Module : mmdb_coormngr <interface>
+// ~~~~~~~~~
+// Project : MacroMolecular Data Base (MMDB)
+// ~~~~~~~~~
+// **** Classes : mmdb::Brick ( space brick )
+// ~~~~~~~~~ mmdb::CoorManager ( MMDB atom coordinate manager )
+//
+// (C) E. Krissinel 2000-2013
+//
+// =================================================================
+//
+
+#ifndef __MMDB_CoorMngr__
+#define __MMDB_CoorMngr__
+
+#include "mmdb_root.h"
+
+namespace mmdb {
+
+ // =========================== Brick ==============================
+
+ // bricking control
+ enum BRICK_STATE {
+ BRICK_ON_1 = 0x00000001,
+ BRICK_ON_2 = 0x00000002,
+ BRICK_READY = 0x00000004
+ };
+
+ DefineClass(Brick);
+ typedef PPBrick * PPPBrick;
+
+ class Brick {
+
+ public :
+ int nAtoms; // number of atoms hit into brick
+ PPAtom atom; // pointers to atoms
+ ivector id; // atom ids (in present realization, these are
+ // indices of atoms from the bricked array)
+
+ Brick ();
+ ~Brick();
+
+ void Clear ();
+ void AddAtom ( PAtom A, int atomid );
+
+ protected :
+ int nAllocAtoms;
+ void InitBrick();
+
+ };
+
+
+ // =========================== MBrick =============================
+
+ // Bricking multiple structures
+
+ DefineClass(MBrick);
+ typedef PPMBrick * PPPMBrick;
+
+ class MBrick {
+
+ public :
+ ivector nAtoms; // number of atoms in the brick
+ PPAtom *atom; // pointers to atoms
+ imatrix id; // atom ids (in present realization, these are
+ // indices of atoms from the bricked array)
+
+ MBrick ( int nStructures );
+ ~MBrick();
+
+ void Clear ();
+ void AddAtom ( PAtom A, int structNo, int atomid );
+
+ protected :
+ ivector nAlloAtoms;
+ int nStruct;
+ void InitMBrick ( int nStructures );
+
+ };
+
+
+
+ // ==================== GenSym ========================
+
+ DefineClass(GenSym);
+ DefineStreamFunctions(GenSym);
+
+ class GenSym : public SymOps {
+
+ friend class CoorManager;
+
+ public :
+
+ GenSym ();
+ GenSym ( io::RPStream Object );
+ ~GenSym();
+
+ void FreeMemory();
+
+ int AddSymOp ( cpstr XYZOperation );
+ // the number of just added operation may be obtained as
+ // Nop = GenSym::GetNofSymOps()-1 .
+
+ int AddRenChain ( int Nop, const ChainID ch1, const ChainID ch2 );
+
+ void Copy ( PSymOps genSym );
+
+ void write ( io::RFile f );
+ void read ( io::RFile f );
+
+ protected :
+
+ PChainID * chID1; // pairs of chains to rename from chID1[n][i]
+ PChainID * chID2; // to chID2[n][i] for each operation n<Nops
+ ivector nChains; // number of chains to rename for each oper-n
+
+ void InitGenSym();
+
+ private :
+ int nOpAlloc; // number of allocated operations
+
+ };
+
+
+ // ========================= Contact =============================
+
+ DefineStructure(Contact);
+
+ struct Contact {
+ int id1,id2;
+ long group;
+ realtype dist;
+ void Copy ( RContact c );
+ void Swap ( RContact c );
+ };
+
+
+ // ======================== MContact =============================
+
+ DefineClass(MContact);
+
+ class MContact : public io::Stream {
+
+ public :
+ int nStruct,contactID;
+ ivector nAtoms;
+ PPAtom * atom;
+ imatrix id;
+
+ MContact ( int nStructures );
+ ~MContact();
+
+ void AddContact ( PAtom A, int structNo, int atomid );
+
+ protected:
+ ivector nAlloc;
+
+ };
+
+ extern void DeleteMContacts ( PPMContact & mcontact, int nContacts );
+
+
+ // ====================== CoorManager =========================
+
+ DefineClass(CoorManager);
+ DefineStreamFunctions(CoorManager);
+
+ // ---- Atom extraction return
+ enum CID_RC {
+ CID_Ok = 0,
+ CID_NoModel = 1,
+ CID_NoChain = 2,
+ CID_NoResidue = 3,
+ CID_NoAtom = 4,
+ CID_WrongPath = 5
+ };
+
+ // ---- generate symmetry mates return codes
+ enum GSM_RC {
+ GSM_Ok = 0,
+ GSM_NoSymOps = 1,
+ GSM_NoTransfMatrices = 2,
+ GSM_NoCell = 3
+ };
+
+ class CoorManager : public Root {
+
+ public :
+
+ int CoorIDCode; // last return from atom extraction procedure
+
+ CoorManager ();
+ CoorManager ( io::RPStream Object );
+ ~CoorManager();
+
+
+ // ----------------------------------------------------------
+
+ int SetDefaultCoorID ( cpstr CID );
+
+
+ // ---------------- Bricking ------------------------------
+
+ void RemoveBricks ();
+ bool areBricks () { return (brick!=NULL); }
+ void MakeBricks ( PPAtom atmvec, int avlen,
+ realtype Margin, realtype BrickSize=6.0 );
+ void GetBrickDimension (
+ int & nxmax, int & nymax, int & nzmax );
+ void GetBrickCoor ( PAtom A, int & nx, int & ny, int & nz );
+ void GetBrickCoor ( realtype x, realtype y, realtype z,
+ int & nx, int & ny, int & nz );
+ PBrick GetBrick ( int nx, int ny, int nz );
+
+ void RemoveMBricks ();
+ bool areMBricks () { return (mbrick!=NULL); }
+ void MakeMBricks ( PPAtom * atmvec, ivector avlen,
+ int nStructures, realtype Margin,
+ realtype BrickSize=6.0 );
+ void GetMBrickDimension (
+ int & nxmax, int & nymax, int & nzmax );
+ void GetMBrickCoor ( PAtom A, int & nx, int & ny, int & nz );
+ void GetMBrickCoor ( realtype x, realtype y, realtype z,
+ int & nx, int & ny, int & nz );
+ PMBrick GetMBrick ( int nx, int ny, int nz );
+
+ // ---------------- Extracting models ---------------------
+
+ int GetNumberOfModels () { return nModels; }
+ int GetFirstModelNum ();
+ PModel GetFirstDefinedModel();
+ PModel GetModel ( int modelNo ); // 1<=modelNo<=nModels
+ PModel GetModel ( cpstr CID );
+ void GetModelTable ( PPModel & modTable,
+ int & NumberOfModels );
+
+ // ---------------- Deleting models -----------------------
+
+ int DeleteModel ( cpstr CID );
+ int DeleteModel ( int modelNo ); // 1<=modelNo<=nOfModels
+
+ // ---------------- Adding/Inserting models ---------------
+
+ int AddModel ( PModel mdl );
+ int InsModel ( PModel mdl, int modelNo );
+ void RotateModels ( int modelNo1, int modelNo2, int rotdir );
+ void SwapModels ( int modelNo1, int modelNo2 );
+
+ // ---------------- Extracting chains ---------------------
+
+ int GetNumberOfChains ( int modelNo );
+ int GetNumberOfChains ( cpstr CID );
+ PChain GetChain ( int modelNo, const ChainID chainID );
+ PChain GetChain ( int modelNo, int chainNo );
+ PChain GetChain ( cpstr CID );
+ void GetChainTable ( int modelNo, PPChain & chainTable,
+ int & NumberOfChains );
+ void GetChainTable ( cpstr CID, PPChain & chainTable,
+ int & NumberOfChains );
+
+ // ----------------- Deleting chains ----------------------
+
+ int DeleteChain ( int modelNo, const ChainID chID );
+ int DeleteChain ( int modelNo, int chainNo );
+ int DeleteAllChains ( int modelNo );
+ int DeleteAllChains ();
+
+ // ------------------ Adding chains -----------------------
+
+ int AddChain ( int modelNo, PChain chain );
+
+ // ---------------- Extracting residues -------------------
+
+ int GetNumberOfResidues ( int modelNo, const ChainID chainID );
+ int GetNumberOfResidues ( int modelNo, int chainNo );
+ int GetNumberOfResidues ( cpstr CID );
+ PResidue GetResidue ( int modelNo, const ChainID chainID,
+ int seqNo, const InsCode insCode );
+ PResidue GetResidue ( int modelNo, int chainNo,
+ int seqNo, const InsCode insCode );
+ PResidue GetResidue ( int modelNo, const ChainID chainID,
+ int resNo );
+ PResidue GetResidue ( int modelNo, int chainNo, int resNo );
+ PResidue GetResidue ( cpstr CID );
+ int GetResidueNo ( int modelNo, const ChainID chainID,
+ int seqNo, const InsCode insCode );
+ int GetResidueNo ( int modelNo, int chainNo,
+ int seqNo, const InsCode insCode );
+ void GetResidueTable ( PPResidue & resTable,
+ int & NumberOfResidues );
+ void GetResidueTable ( int modelNo, const ChainID chainID,
+ PPResidue & resTable,
+ int & NumberOfResidues );
+ void GetResidueTable ( int modelNo, int chainNo,
+ PPResidue & resTable,
+ int & NumberOfResidues );
+ void GetResidueTable ( cpstr CID, PPResidue & resTable,
+ int & NumberOfResidues );
+
+
+ // ----------------- Deleting residues -----------------------
+
+ int DeleteResidue ( int modelNo, const ChainID chainID,
+ int seqNo, const InsCode insCode );
+ int DeleteResidue ( int modelNo, const ChainID chainID,
+ int resNo );
+ int DeleteResidue ( int modelNo, int chainNo,
+ int seqNo, const InsCode insCode );
+ int DeleteResidue ( int modelNo, int chainNo, int resNo );
+ int DeleteAllResidues ( int modelNo, const ChainID chainID );
+ int DeleteAllResidues ( int modelNo, int chainNo );
+ int DeleteAllResidues ( int modelNo );
+ int DeleteAllResidues ();
+ int DeleteSolvent ();
+
+ // ------------------- Adding residues -----------------------
+
+ int AddResidue ( int modelNo, const ChainID chainID,
+ PResidue res );
+ int AddResidue ( int modelNo, int chainNo, PResidue res );
+
+ // -------------------- Extracting atoms ----------------------
+
+ int GetNumberOfAtoms () { return nAtoms; }
+ int GetNumberOfAtoms ( int modelNo, const ChainID chainID,
+ int seqNo, const InsCode insCode );
+ int GetNumberOfAtoms ( int modelNo, int chainNo,
+ int seqNo, const InsCode insCode );
+ int GetNumberOfAtoms ( int modelNo, const ChainID chainID,
+ int resNo );
+ int GetNumberOfAtoms ( int modelNo, int chainNo, int resNo );
+ int GetNumberOfAtoms ( cpstr CID );
+
+ PAtom GetAtom (
+ int modelNo, // model serial number 1...
+ const ChainID chID, // chain ID
+ int seqNo, // residue sequence number
+ const InsCode insCode, // residue insertion code
+ const AtomName aname, // atom name
+ const Element elmnt, // chemical element code or '*'
+ const AltLoc aloc // alternate location indicator
+ );
+
+ PAtom GetAtom (
+ int modelNo, // model serial number 1...
+ const ChainID chID, // chain ID
+ int seqNo, // residue sequence number
+ const InsCode insCode, // residue insertion code
+ int atomNo // atom number 0..
+ );
+
+ PAtom GetAtom (
+ int modelNo, // model serial number 1...
+ const ChainID chID, // chain ID
+ int resNo, // residue number 0..
+ const AtomName aname, // atom name
+ const Element elmnt, // chemical element code or '*'
+ const AltLoc aloc // alternate location indicator
+ );
+
+ PAtom GetAtom (
+ int modelNo, // model serial number 1...
+ const ChainID chID, // chain ID
+ int resNo, // residue number 0..
+ int atomNo // atom number 0..
+ );
+
+ PAtom GetAtom (
+ int modelNo, // model serial number 1...
+ int chNo, // chain number 0..
+ int seqNo, // residue sequence number
+ const InsCode insCode, // residue insertion code
+ const AtomName aname, // atom name
+ const Element elmnt, // chemical element code or '*'
+ const AltLoc aloc // alternate location indicator
+ );
+
+ PAtom GetAtom (
+ int modelNo, // model serial number 1...
+ int chNo, // chain number 0...
+ int seqNo, // residue sequence number
+ const InsCode insCode, // residue insertion code
+ int atomNo // atom number 0...
+ );
+
+ PAtom GetAtom (
+ int modelNo, // model serial number 1...
+ int chNo, // chain number 0...
+ int resNo, // residue number 0...
+ const AtomName aname, // atom name
+ const Element elmnt, // chemical element code or '*'
+ const AltLoc aloc // alternate location indicator
+ );
+
+ PAtom GetAtom (
+ int modelNo, // model serial number 1...
+ int chNo, // chain number 0...
+ int resNo, // residue number 0...
+ int atomNo // atom number 0...
+ );
+
+
+ // GetAtom(CID) returns atom answering to the following
+ // CID pattern:
+ // /mdl/chn/seq(res).i/atm[elm]:a
+ // where
+ // mdl - model number (mandatory); at least model #1 is always
+ // present
+ // chn - chain identifier ( mandatory)
+ // seq - residue sequence number (mandatory)
+ // (res) - residue name in round brackets (may be omitted)
+ // .i - insert code after a dot; if '.i' or 'i' is missing
+ // then residue without an insertion code is looked
+ // for
+ // atm - atom name (mandatory)
+ // [elm] - chemical element code in square brackets; it may
+ // be omitted but could be helpful for e.g.
+ // distinguishing C_alpha and CA
+ // :a - alternate location indicator after colon; if
+ // ':a' or 'a' is missing then an atom without
+ // alternate location indicator is looked for.
+ // All spaces are ignored, all identifiers should be in capital
+ // letters (comparisons are case-sensitive).
+ PAtom GetAtom ( cpstr CID );
+
+
+ void GetAtomTable ( PPAtom & atomTable, int & NumberOfAtoms );
+ void GetAtomTable ( int modelNo, const ChainID chainID,
+ int seqNo, const InsCode insCode,
+ PPAtom & atomTable, int & NumberOfAtoms );
+ void GetAtomTable ( int modelNo, int chainNo,
+ int seqNo, const InsCode insCode,
+ PPAtom & atomTable, int & NumberOfAtoms );
+ void GetAtomTable ( int modelNo, const ChainID chainID, int resNo,
+ PPAtom & atomTable, int & NumberOfAtoms );
+ void GetAtomTable ( int modelNo, int chainNo, int resNo,
+ PPAtom & atomTable, int & NumberOfAtoms );
+ void GetAtomTable ( cpstr CID, PPAtom & atomTable,
+ int & NumberOfAtoms );
+
+
+ // GetAtomTable1(..) returns atom table without TER atoms and
+ // without NULL atom pointers. NumberOfAtoms returns the actual
+ // number of atom pointers in atomTable.
+ // atomTable is allocated within the function. If it was
+ // not set to NULL before calling the function, the function will
+ // attempt to deallocate it first.
+ // The application is responsible for deleting atomTable,
+ // however it must not touch atom pointers, i.e. use simply
+ // "delete atomTable;". Never pass atomTable from GetAtomTable(..)
+ // into this function, unless you set it to NULL before doing that.
+ void GetAtomTable1 ( PPAtom & atomTable, int & NumberOfAtoms );
+ void GetAtomTable1 ( int modelNo, const ChainID chainID,
+ int seqNo, const InsCode insCode,
+ PPAtom & atomTable, int & NumberOfAtoms );
+ void GetAtomTable1 ( int modelNo, int chainNo,
+ int seqNo, const InsCode insCode,
+ PPAtom & atomTable, int & NumberOfAtoms );
+ void GetAtomTable1 ( int modelNo, const ChainID chainID, int resNo,
+ PPAtom & atomTable, int & NumberOfAtoms );
+ void GetAtomTable1 ( int modelNo, int chainNo, int resNo,
+ PPAtom & atomTable, int & NumberOfAtoms );
+ void GetAtomTable1 ( cpstr CID, PPAtom & atomTable,
+ int & NumberOfAtoms );
+
+
+ // -------------------- Deleting atoms -----------------------
+
+ int DeleteAtom ( int modelNo,
+ const ChainID chID,
+ int seqNo,
+ const InsCode insCode,
+ const AtomName aname,
+ const Element elmnt,
+ const AltLoc aloc );
+ int DeleteAtom ( int modelNo,
+ const ChainID chID,
+ int seqNo,
+ const InsCode insCode,
+ int atomNo );
+ int DeleteAtom ( int modelNo,
+ const ChainID chID,
+ int resNo,
+ const AtomName aname,
+ const Element elmnt,
+ const AltLoc aloc );
+ int DeleteAtom ( int modelNo, const ChainID chID,
+ int resNo, int atomNo );
+ int DeleteAtom ( int modelNo, int chNo, int seqNo,
+ const InsCode insCode,
+ const AtomName aname,
+ const Element elmnt,
+ const AltLoc aloc );
+ int DeleteAtom ( int modelNo, int chNo, int seqNo,
+ const InsCode insCode, int atomNo );
+ int DeleteAtom ( int modelNo, int chNo, int resNo,
+ const AtomName aname,
+ const Element elmnt,
+ const AltLoc aloc );
+ int DeleteAtom ( int modelNo, int chNo, int resNo, int atomNo );
+
+ int DeleteAllAtoms ( int modelNo, const ChainID chID,
+ int seqNo, const InsCode insCode );
+ int DeleteAllAtoms ( int modelNo, const ChainID chID, int resNo );
+ int DeleteAllAtoms ( int modelNo, const ChainID chID );
+ int DeleteAllAtoms ( int modelNo, int chNo, int seqNo,
+ const InsCode insCode );
+ int DeleteAllAtoms ( int modelNo, int chNo, int resNo );
+ int DeleteAllAtoms ( int modelNo, int chNo );
+ int DeleteAllAtoms ( int modelNo );
+ int DeleteAllAtoms ();
+
+ // This function leaves only alternative location with maximal
+ // occupancy, if those are equal or unspecified, the one with
+ // "least" alternative location indicator.
+ // The function returns the number of deleted atoms and optimizes
+ // the atom index.
+ int DeleteAltLocs ();
+
+
+ // --------------------- Adding atoms ------------------------
+
+ int AddAtom ( int modelNo, const ChainID chID,
+ int seqNo, const InsCode insCode, PAtom atom );
+ int AddAtom ( int modelNo, const ChainID chID, int resNo,
+ PAtom atom );
+ int AddAtom ( int modelNo, int chNo, int seqNo,
+ const InsCode insCode, PAtom atom );
+ int AddAtom ( int modelNo, int chNo, int resNo, PAtom atom );
+
+
+ // -------------------- Transformations -----------------------
+
+ int GenerateSymMates ( PGenSym genSym=NULL );
+ // 1: no Sym operations,
+ // 2: no fract/orth matrices
+ // 3: no cell parameters
+ // 0: Ok
+
+ void ApplyTransform ( mat44 & TMatrix ); // simply transforms all
+ // coordinates by multiplying
+ // with matrix TMatrix
+
+ int BringToUnitCell(); // brings all chains into 0th unit cell
+
+ // Frac2Orth(..) and Orth2Frac(..) transform between fractional
+ // and orthogonal coordinates, if areMatrices() returns true.
+ // If the transformation matrices were not set, the functions just
+ // copy the coordinates. Returns true if the transformation was
+ // done; False return means that transformation matrices were not
+ // calculated
+ bool Frac2Orth (
+ realtype xfrac, realtype yfrac, realtype zfrac,
+ realtype & xorth, realtype & yorth, realtype & zorth );
+ bool Orth2Frac (
+ realtype xorth, realtype yorth, realtype zorth,
+ realtype & xfrac, realtype & yfrac, realtype & zfrac );
+
+
+ // Below, F and T are transformation matrices in fractional and
+ // orthogonal coordinates, respectively.
+ bool Frac2Orth ( mat44 & F, mat44 & T );
+ bool Orth2Frac ( mat44 & T, mat44 & F );
+
+ // ==================== Seeking contacts ======================
+
+ void SeekContacts (
+ PPAtom AIndex, // index of atoms [0..ilen-1]
+ int ilen, // length of index
+ int atomNum, // number of 1st contact atom
+ // in the index. All other atoms
+ // are checked for contact with
+ // 1st atom
+ realtype dist1, // minimal contact distance
+ realtype dist2, // maximal contact distance
+ int seqDist, // the sequence distance to neglect.
+ // If seqDist==0, all atoms are
+ // checked for contact. If
+ // seqDist==1, the atoms belonging
+ // to the same residue as atom
+ // AIndex[atomNum], are neglected.
+ // If seqDist>1, all atoms belonging
+ // to residues closer than
+ // +/-(seqDist-1) around that of
+ // atom AIndex[atomNum], are
+ // neglected. If chain is broken
+ // (has a gap) on section
+ // [-(seqDist-1)..seqDist-1], the
+ // section of neglection is
+ // shortened to that gap.
+ RPContact contact, // indices of contacting atoms
+ // [0..ncontacts-1]. contact[i].id1
+ // is set to atomNum and
+ // contact[i].id2 is set to the
+ // index of 2nd contacting atom
+ // in vector AIndex
+ int & ncontacts, // number of contacts found. If
+ // ncontacts>0 on input, it is
+ // assumed that new contacts that
+ // newly found contacts should be
+ // appended to those already
+ // existing
+ int maxlen=0, // if <=0, then vector contact is
+ // allocated dynamically. If
+ // contact!=NULL, then it is
+ // appended with new contacts.
+ // The application is responsible
+ // for deallocation of contact
+ // after use.
+ // If maxlen>0 then vector contact
+ // is prohibited of dynamical
+ // allocation/deallocation. In this
+ // case, not more than maxlen
+ // contacts will be returned.
+ long group=0 // a contact group ID, which will be
+ // simply stored in contact[i].group
+ // fields. This ID may be useful
+ // if contacts are obtained in
+ // multiple calls of the function
+ );
+
+ void SeekContacts (
+ PAtom A, // 1st atom in contact
+ PPAtom AIndex, // index of atoms [0..ilen-1] to
+ // check for contact with 1st atom
+ int ilen, // length of index
+ realtype dist1, // minimal contact distance
+ realtype dist2, // maximal contact distance
+ int seqDist, // the sequence distance to neglect.
+ // If seqDist==0, all atoms are
+ // checked for contact. If
+ // seqDist==1, the atoms belonging
+ // to the same residue as atom
+ // A, are neglected. If seqDist>1,
+ // all atoms belonging to residues
+ // closer than +/-(seqDist-1) around
+ // that of atom A, are neglected. If
+ // chain is broken (has a gap) on
+ // section
+ // [-(seqDist-1)..seqDist-1], the
+ // section of neglection is
+ // shortened to that gap.
+ RPContact contact, // indices of contacting atoms
+ // [0..ncontacts-1]. contact[i].id1
+ // is set to -1, and contact[i].id2
+ // is set to the index of 2nd
+ // contacting atom in vector AIndex
+ int & ncontacts, // number of contacts found. If
+ // ncontacts>0 on input, it is
+ // assumed that new contacts that
+ // newly found contacts should be
+ // appended those already existing
+ int maxlen=0, // if <=0, then vector contact is
+ // allocated dynamically. If
+ // contact!=NULL, then it is
+ // appended with new contacts.
+ // The application is responsible
+ // for deallocation of contact
+ // after use.
+ // If maxlen>0 then vector contact
+ // is prohibited of dynamical
+ // allocation/deallocation. In this
+ // case, not more than maxlen
+ // contacts will be returned.
+ long group=0 // a contact group ID, which will be
+ // simply stored in contact[i].group
+ // fields. This ID may be useful
+ // if contacts are obtained in
+ // multiple calls of the function
+ );
+
+ void SeekContacts (
+ PPAtom AIndex1, // 1st atom index [0..ilen1-1]
+ int ilen1, // length of 1st index
+ PPAtom AIndex2, // 2nd atom index [0..ilen2-1] to
+ // check for contact with 1st index
+ int ilen2, // length of 2nd index
+ realtype dist1, // minimal contact distance
+ realtype dist2, // maximal contact distance
+ int seqDist, // the sequence distance to
+ // neglect.
+ // If seqDist==0, all atoms are
+ // checked for contact.
+ // If seqDist==1, the atoms
+ // belonging to the same residue
+ // are neglected.
+ // If seqDist>1, all atoms
+ // belonging to residues closer than
+ // +/-(seqDist-1) to each other,
+ // are neglected. If chain is broken
+ // (has a gap) on section
+ // [-(seqDist-1)..seqDist-1], the
+ // section of neglection is
+ // shortened to that gap.
+ RPContact contact, // indices of contacting atoms
+ // [0..ncontacts-1]. contact[i].id1
+ // contains number of atom from 1st
+ // index, and contact[i].id2
+ // contains number of atom from 2nd
+ // index, contacting with the former
+ // one
+ int & ncontacts, // number of contacts found. If
+ // ncontacts>0 on input, it is
+ // assumed that newly found
+ // contacts should be appended to
+ // those already existing
+ int maxlen=0, // if <=0, then vector contact is
+ // allocated dynamically. If
+ // contact!=NULL, then it is
+ // appended with new contacts.
+ // The application is responsible
+ // for deallocation of contact
+ // after use.
+ // If maxlen>0 then vector contact
+ // is prohibited of dynamical
+ // allocation/deallocation. In this
+ // case, not more than maxlen
+ // contacts will be returned.
+ mat44 * TMatrix=NULL, // transformation matrix for 2nd
+ // set of atoms (AIndex2)
+ long group=0, // a contact group ID, which will
+ // be stored in contact[i].group
+ // fields. This ID may be useful
+ // if contacts are obtained in
+ // multiple calls of the function
+ int bricking=0, // bricking control; may be a
+ // combination of BRICK_ON_1 or
+ // BRICK_ON_2 with BRICK_READY
+ bool doSqrt=true // if False, then Contact contains
+ // square distances
+ );
+
+ // Simplified optimized for speed version:
+ // - no NULL pointers and Ters in AIndex1 and AIndex2
+ // - no checks for identity atoms in AIndex1 and AIndex2
+ // - contact must be pre-allocated with at least ilen1*ilen2
+ // elements
+ // - contact returns square distances
+ // - ncontacts is always reset
+ void SeekContacts (
+ PPAtom AIndex1, // 1st atom index [0..ilen1-1]
+ int ilen1, // length of 1st index
+ PPAtom AIndex2, // 2nd atom index [0..ilen2-1] to
+ // check for contact with 1st index
+ int ilen2, // length of 2nd index
+ realtype contDist, // maximal contact distance
+ PContact contact, // indices of contacting atoms
+ // [0..ncontacts-1]. contact[i].id1
+ // contains number of atom from 1st
+ // index, and contact[i].id2
+ // contains number of atom from 2nd
+ // index, contacting with the former
+ // one. Must be pre-allocated
+ int & ncontacts, // number of contacts found
+ int bricking=0 // bricking control; may be a
+ // combination of BRICK_ON_1 or
+ // BRICK_ON_2 with BRICK_READY
+ );
+
+ void SeekContacts (
+ PPAtom AIndex1, // 1st atom index [0..ilen1-1]
+ int ilen1, // length of 1st index
+ PPAtom * AIndex2, // indexes of atoms to be checked
+ // for contact with each atom from
+ // Aindex1; dimension
+ // [0..nStructures-1][0..ilen2[i]-1]
+ ivector ilen2, // lengths of indexes AIndex2
+ int nStructures, // number of indexes AIndex2
+ realtype dist1, // minimal contact distance
+ realtype dist2, // maximal contact distance
+ PPMContact & contact, // resulting contacts, one structure
+ // per each position in AIndex1. If
+ // AIndex1[i] is NULL, contact[i] is
+ // also NULL. "contact" is always
+ // allocated, no re-use or
+ // re-allocation is attempted.
+ int bricking=0 // bricking control; may be
+ // BRICK_READY if AIndex2 does not
+ // change
+ );
+
+ protected :
+
+ // bricks
+ realtype brick_size, xbrick_0,ybrick_0,zbrick_0;
+ int nbrick_x,nbrick_y,nbrick_z;
+ PPPBrick * brick;
+
+ realtype mbrick_size, xmbrick_0,ymbrick_0,zmbrick_0;
+ int nmbrick_x,nmbrick_y,nmbrick_z;
+ PPPMBrick * mbrick;
+
+ // --------------- Stream I/O -----------------------------
+ void write ( io::RFile f );
+ void read ( io::RFile f );
+
+ void InitMMDBCoorManager();
+
+ void ApplySymTransform ( int SymMatrixNo, PGenSym genSym=NULL );
+
+ void ResetManager ();
+
+ void FindSeqSection ( PAtom atom, int seqDist,
+ int & seq1, int & seq2 );
+ bool iContact ( PAtom a1, PAtom a2,
+ int seq1, int seq2,
+ realtype dd, realtype d12,
+ realtype d22, realtype & d2 );
+ bool iContact ( realtype x, realtype y,
+ realtype z, PAtom a2,
+ realtype dd, realtype d12,
+ realtype d22, realtype & d2 );
+
+ };
+
+
+
+ // ===================================================================
+
+
+
+ // GetEulerRotMatrix(..) calculates the Euler rotation matrix
+ // for rotation:
+ // 1) about z-axis by angle alpha
+ // 2) about new y-axis by angle beta
+ // 3) about new z-axis by angle gamma
+ extern void GetEulerRotMatrix ( mat33 & erm, realtype alpha,
+ realtype beta, realtype gamma );
+
+ // GetEulerTMatrix(..) calculates the Euler rotation-translation
+ // matrix for rotation:
+ // 1) about z-axis by angle alpha
+ // 2) about new y-axis by angle beta
+ // 3) about new z-axis by angle gamma
+ // Point (x0,y0,z0) is the center of rotation.
+ extern void GetEulerTMatrix ( mat44 & erm, realtype alpha,
+ realtype beta, realtype gamma,
+ realtype x0, realtype y0, realtype z0 );
+
+ // Euler rotation: 1) about z-axis by angle alpha
+ // 2) about new y-axis by angle beta
+ // 3) about new z-axis by angle gamma
+ // Point (x0,y0,z0) is the center of rotation.
+ extern void EulerRotation ( PPAtom A, int nA,
+ realtype alpha, realtype beta, realtype gamma,
+ realtype x0, realtype y0, realtype z0 );
+
+ // GetVecRotMatrix(..) calculates the rotation matrix for
+ // rotation by angle alpha about arbitrary vector directed
+ // as (vx,vy,vz) = (vx2-vx1,vy2-vy1,vz2-vz1).
+ extern void GetVecRotMatrix ( mat33 & vrm, realtype alpha,
+ realtype vx, realtype vy, realtype vz );
+
+
+ // Given the rotation matrix vrm, GetRotParameters(..)
+ // returns the rotation angle alpha and the normalized
+ // rotation axis vector (vx,vy,vz).
+ // The rotation angle and vector are determined up to
+ // their sign (however correlated, so that being substituted
+ // into GetVecRotMatrix(..) they yield the same rotation
+ // matrix).
+ // The function does not check for vrm to be a valid
+ // rotation matrix.
+ extern void GetRotParameters ( mat33 & vrm, realtype & alpha,
+ realtype & vx, realtype & vy, realtype & vz );
+
+
+ // GetVecTMatrix(..) calculates the rotation-translation matrix
+ // for rotation by angle alpha about arbitrary vector directed as
+ // (vx,vy,vz) = (vx2-vx1,vy2-vy1,vz2-vz1). Point (x0,y0,z0) is
+ // the center of rotation -- actually a point belonging to the
+ // rotation axis.
+ extern void GetVecTMatrix ( mat44 & vrm, realtype alpha,
+ realtype vx, realtype vy, realtype vz,
+ realtype x0, realtype y0, realtype z0 );
+
+ // Vector rotation is rotation by angle alpha about arbitrary
+ // vector directed as (vx,vy,vz) = (vx2-vx1,vy2-vy1,vz2-vz1).
+ // Point (x0,y0,z0) is the center of rotation -- actually
+ // a point belonging to the rotation axis.
+ extern void VectorRotation ( PPAtom A, int nA, realtype alpha,
+ realtype vx, realtype vy, realtype vz,
+ realtype x0, realtype y0, realtype z0 );
+
+ extern void GetMassCenter ( PPAtom A, int nA,
+ realtype & xmc, realtype & ymc, realtype & zmc );
+
+
+ enum SPOSEAT_RC {
+ SPOSEAT_Ok = 0,
+ SPOSEAT_NoAtoms = 1,
+ SPOSEAT_SVD_Fail = 2
+ };
+
+ // Given two sets of atoms, A1 and A2, SuperposeAtoms(...) calculates
+ // the rotational-translational matrix T such that |T*A1 - A2| is
+ // minimal in least-square terms.
+ // If vector C is not given (default), all nA atoms of set A1 are
+ // considered as corresponding to nA first atoms of set A2,
+ // A1[i] <-> A2[i], 0<=i<nA .
+ // If vector C is given, then the correspondence of atoms is
+ // established as A1[i] <-> A2[C[i]] only for those i that C[i]>=0.
+ // The default option (C==NULL) is thus identical to C[i]==i, 0<=i<nA.
+ // Upon normal completion, the procedure returns SPOSEAT_Ok.
+
+ extern int SuperposeAtoms ( mat44 & T, PPAtom A1, int nA, PPAtom A2,
+ ivector C=NULL );
+
+ enum CNSORT_DIR {
+ CNSORT_OFF = 0,
+ CNSORT_1INC = 1,
+ CNSORT_1DEC = 2,
+ CNSORT_2INC = 3,
+ CNSORT_2DEC = 4,
+ CNSORT_DINC = 5,
+ CNSORT_DDEC = 6
+ };
+
+ extern void SortContacts ( PContact contact, int ncontacts,
+ CNSORT_DIR sortmode );
+
+
+ extern const realtype NO_TORSION;
+
+ extern realtype getPhi ( PPAtom A ); // A[0] - A[3] used
+ extern realtype getPsi ( PPAtom A ); // A[0] - A[2] used
+
+} // namespace mmdb
+
+#endif
+