pdb.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors:     Carlos Oscar Sanchez Sorzano (coss@cnb.csic.es)
 *
 * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * 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 General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
 * 02111-1307  USA
 *
 *  All comments concerning this program package may be sent to the
 *  e-mail address 'xmipp@cnb.csic.es'
 ***************************************************************************/
/*****************************************************************************/
/* INTERACTION WITH PDBs                                                     */
/*****************************************************************************/

#ifndef _XMIPP_PDB_HH
#define _XMIPP_PDB_HH

#include <vector>
#include "cif++.hpp"
#include "core/xmipp_error.h"

template<typename T>
class Matrix1D;
template<typename T>
class Matrix2D;
template<typename T>
class MultidimArray;
class FileName;
class Projection;
class Histogram1D;

int atomCharge(const std::string &atom);

double atomRadius(const std::string &atom);

double atomCovalentRadius(const std::string &atom);

void computePDBgeometry(const std::string &fnPDB,
                        Matrix1D<double> &centerOfMass,
                        Matrix1D<double> &limit0, Matrix1D<double> &limitF,
                        const std::string &intensityColumn);

void applyGeometryToPDBFile(const std::string &fn_in, const std::string &fn_out,
                   const Matrix2D<double> &A, bool centerPDB=true,
                   const std::string &intensityColumn="occupancy");

struct pdbInfo
{
    std::vector<double> x;
    std::vector<double> y;
    std::vector<double> z;
    std::vector<double> b;
    std::vector<std::string> chain;
    std::vector<int> residue;
    std::vector<double> atomCovRad;
};

void analyzePDBAtoms(const FileName &fn_pdb, const std::string &typeOfAtom, int &numberOfAtoms, pdbInfo &at_pos);


class Atom
{
public:
    char atomType;

    double x;

    double y;

    double z;
};

class PDBPhantom
{
public:
    std::vector<Atom> atomList;

    // Whole data block
    cif::datablock dataBlock;

    void addAtom(const Atom &atom)
    {
        atomList.push_back(atom);
    }

    const Atom& getAtom(int i) const
    {
        return atomList[i];
    }

    size_t getNumberOfAtoms() const
    {
        return atomList.size();
    }

    void read(const FileName &fnPDB);

    void shift(double x, double y, double z);

    void produceSideInfo();
};

class RichAtom
{
public:
    std::string record;

    int serial;

    double x;

    double y;

    double z;

    std::string name;

    std::string altloc;

    std::string resname;

    char chainid;

    int resseq;

    std::string icode;

    double occupancy;

    double bfactor;

    std::string atomType;

    std::string charge;

    /* PDB Specific values */
    std::string segment;

    /* CIF Specific values */
    // Alternative id
    std::string altId;

    // Sequence id
    int seqId;

    // Author sequence id
    int authSeqId;

    // Author chain name
    std::string authCompId;

    // Author chain location
    std::string authAsymId;

    // Author atom name
    std::string authAtomId;

    // PDB model number
    int pdbNum;
};

class PDBRichPhantom
{
public:
    std::vector<std::string> remarks;

    std::vector<RichAtom> atomList;
    std::vector<double> intensities;

    // Whole data block
    cif::datablock dataBlock;

    void addAtom(const RichAtom &atom)
    {
        atomList.push_back(atom);
    }

    size_t getNumberOfAtoms() const
    {
        return atomList.size();
    }

    void read(const FileName &fnPDB, const bool pseudoatoms = false, const double threshold = 0.0);

    void write(const FileName &fnPDB, const bool renumber = false);

};

void atomDescriptors(const std::string &atom, Matrix1D<double> &descriptors);

double electronFormFactorFourier(double f,
                                 const Matrix1D<double> &descriptors);

double electronFormFactorRealSpace(double r,
                                   const Matrix1D<double> &descriptors);

void atomRadialProfile(int M, double T, const std::string &atom,
                       Matrix1D<double> &profile);

void atomProjectionRadialProfile(int M,
                                 const Matrix1D<double> &profileCoefficients,
                                 Matrix1D<double> &projectionProfile);

class AtomInterpolator
{
public:
    // Vector of radial volume profiles
    std::vector< MultidimArray<double> > volumeProfileCoefficients;
    // Vector of radial projection profiles
    std::vector< MultidimArray<double> > projectionProfileCoefficients;
    // Vector of atom radii
    std::vector<double> radii;
    // Downsampling factor
    int M;
    // Fine sampling rate
    double highTs;

    void setup(int m, double hights, bool computeProjection=false);

    void addAtom(const std::string &atomType, bool computeProjection=false);

    int getAtomIndex(char atom) const
    {
        int idx=-1;
        switch (atom)
        {
        case 'H':
            idx=0;
            break;
        case 'C':
            idx=1;
            break;
        case 'N':
            idx=2;
            break;
        case 'O':
            idx=3;
            break;
        case 'P':
            idx=4;
            break;
        case 'S':
            idx=5;
            break;
        case 'E': // Iron Fe
            idx=6;
            break;
        case 'K':
            idx=7;
            break;
        case 'F':
             idx=8;
             break;
        case 'G': // Magnesium Mg
             idx=9;
             break;
        case 'L': // Chlorine Cl
             idx=10;
             break;
        case 'A': // Calcium Ca
             idx=11;
             break;
        default:
            REPORT_ERROR(ERR_VALUE_INCORRECT,(std::string)
                         "AtomInterpolator::getAtomIndex: Atom "+atom+" unknown");
        }
        return idx;
    }

    double atomRadius(char atom) const
    {
        return radii[getAtomIndex(atom)];
    }

    double volumeAtDistance(char atom, double r) const;

    double projectionAtDistance(char atom, double r) const;
};

void projectPDB(const PDBPhantom &phantomPDB,
                const AtomInterpolator &interpolator, Projection &proj,
                int Ydim, int Xdim, double rot, double tilt, double psi);

void distanceHistogramPDB(const PDBPhantom &phantomPDB, size_t Nnearest, double maxDistance, int Nbins, Histogram1D &hist);
#endif

const char*
hy36encode(unsigned width, int value, char* result);

const char*
hy36decode(unsigned width, const char* s, unsigned s_size, int* result);

void
hy36encodeSafe(unsigned width, int value, char* result);

void
hy36decodeSafe(unsigned width, const char* s, unsigned s_size, int* result);