PDB

Collaboration diagram for PDB:

Classes
struct	pdbInfo
class	Atom
class	PDBPhantom
class	RichAtom
class	PDBRichPhantom
class	AtomInterpolator

Functions
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")
void	analyzePDBAtoms (const FileName &fn_pdb, const std::string &typeOfAtom, int &numberOfAtoms, pdbInfo &at_pos)
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)
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)

Detailed Description

Function Documentation

analyzePDBAtoms()

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

ANALYZEPDBDATA takes as input a filename of a pdb file (atomic model) and selects only the typeOfAtom, (for instance the C-alpha atoms) storing the atom positions, b- factor, the residue of each atom, and the covalent radiues in a struct vector at_pos. Also the number of atoms is kept.

Definition at line 64 of file pdb.cpp.

{
    //Open the pdb file
    std::ifstream f2parse;
    f2parse.open(fn_pdb.c_str());

    numberOfAtoms = 0;

    // Initializing and reading pdb file
    PDBRichPhantom pdbFile;

    // Read centered pdb
    pdbFile.read(fn_pdb.c_str());

    // For each atom, store necessary info if type matches
    for (auto& atom : pdbFile.atomList) {
        if (atom.name == typeOfAtom) {
            numberOfAtoms++;

            // Storing coordinates
            at_pos.x.push_back(atom.x);
            at_pos.y.push_back(atom.y);
            at_pos.z.push_back(atom.z);
            at_pos.chain.push_back(std::string(1, atom.chainid));

            // Residue Number
            at_pos.residue.push_back(atom.resseq);

            // Getting the bfactor = 8pi^2*u
            at_pos.b.push_back(atom.bfactor); //sqrt(atom.bfactor/(8*PI*PI));

            // Covalent radius of the atom
            at_pos.atomCovRad.push_back(atomCovalentRadius(atom.name));
        }
    }
}

applyGeometryToPDBFile()

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

Apply geometry transformation to an input PDB. The result is written in the output PDB. Set centerPDB if you want to compute the center of mass first and apply the transformation after centering the PDB.

Definition at line 294 of file pdb.cpp.

{
    Matrix1D<double> centerOfMass;
    Matrix1D<double> limit0;
    Matrix1D<double> limitF;
    if (centerPDB)
    {
        computePDBgeometry(fn_in, centerOfMass,limit0, limitF,
                           intensityColumn);
        limit0 -= centerOfMass;
        limitF -= centerOfMass;
    }

    // Open files
    std::ifstream fh_in;
    fh_in.open(fn_in.c_str());
    if (!fh_in)
        REPORT_ERROR(ERR_IO_NOTEXIST, fn_in);
    std::ofstream fh_out;
    fh_out.open(fn_out.c_str());
    if (!fh_out)
        REPORT_ERROR(ERR_IO_NOWRITE, fn_out);

    // Process all lines of the file
    while (!fh_in.eof())
    {
        // Read an ATOM line
        std::string line;
        getline(fh_in, line);
        if (line == "")
        {
            fh_out << "\n";
            continue;
        }
        std::string kind = line.substr(0,4);
        if (kind != "ATOM" && kind != "HETA")
        {
            fh_out << line << std::endl;
            continue;
        }

        // Extract atom type and position
        // Typical line:
        // ATOM    909  CA  ALA A 161      58.775  31.984 111.803  1.00 34.78
        double x = textToFloat(line.substr(30,8));
        double y = textToFloat(line.substr(38,8));
        double z = textToFloat(line.substr(46,8));

        Matrix1D<double> v(4);
        if (centerPDB)
        {
            VECTOR_R3(v,x-XX(centerOfMass),
                      y-YY(centerOfMass),z-ZZ(centerOfMass));
        }
        else
        {
            VECTOR_R3(v,x,y,z);
        }
        v(3)=1;
        v=A*v;

        char aux[15];
        sprintf(aux,"%8.3f",XX(v));
        line.replace(30,8,aux);
        sprintf(aux,"%8.3f",YY(v));
        line.replace(38,8,aux);
        sprintf(aux,"%8.3f",ZZ(v));
        line.replace(46,8,aux);

        fh_out << line << std::endl;
    }

    // Close files
    fh_in.close();
    fh_out.close();
}

atomCharge()

int atomCharge

(

const std::string &

atom

)

Returns the charge of an atom. Returns 0 if the atom is not within the short list (H, C, N, O, S, P, Fe) of valid atoms.

Definition at line 122 of file pdb.cpp.

{
    switch (atom[0])
    {
    case 'H':
        return 1;
        break;
    case 'C':
        return 6;
        break;
    case 'N':
        return 7;
        break;
    case 'O':
        return 8;
        break;
    case 'P':
        return 15;
        break;
    case 'S':
        return 16;
        break;
    case 'E': // Iron Fe
        return 26;
        break;
    case 'K':
        return 19;
        break;
    case 'F':
        return 9;
        break;
    case 'G': // Magnesium Mg
        return 12;
        break;
    case 'L': // Chlorine Cl
        return 17;
        break;
    case 'A': // Calcium Ca
        return 20;
        break;
    default:
        return 0;
    }
}

atomCovalentRadius()

double atomCovalentRadius

(

const std::string &

atom

)

Returns the covalent radius of an atom. Returns 0 if the atom is not within the short list (H, C, N, O, S, P, Fe) of valid atoms. The radius data is taken from http://www.webelements.com as the empirical radius.

Definition at line 214 of file pdb.cpp.

{
    switch (atom[0])
    {
    case 'H':
        return 0.38;
    case 'C':
        return 0.77;
    case 'N':
        return 0.75;
    case 'O':
        return 0.73;
    case 'P':
        return 1.06;
    case 'S':
        return 1.02;
    case 'F': // Iron
        return 1.25;
    default:
        return 0;
    }
}

atomDescriptors()

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

Description of the electron scattering factors. The returned descriptor is descriptor(0)=Z (number of electrons of the atom), descriptor(1-5)=a1-5, descriptor(6-10)=b1-5. The electron scattering factor at a frequency f (Angstroms^-1) is computed as f_el(f)=sum_i(ai exp(-bi*x^2)). Use the function electronFormFactorFourier or

See Peng, Ren, Dudarev, Whelan. Robust parameterization of elastic and absorptive electron atomic scattering factors. Acta Cryst. A52: 257-276 (1996). Table 3 and equation 3.

Definition at line 883 of file pdb.cpp.

{
    descriptors.initZeros(11);
    if (atom=="H")
    {
        descriptors( 0)= 1;     // Z
        descriptors( 1)= 0.0088; // a1
        descriptors( 2)= 0.0449; // a2
        descriptors( 3)= 0.1481; // a3
        descriptors( 4)= 0.2356; // a4
        descriptors( 5)= 0.0914; // a5
        descriptors( 6)= 0.1152; // b1
        descriptors( 7)= 1.0867; // b2
        descriptors( 8)= 4.9755; // b3
        descriptors( 9)=16.5591; // b4
        descriptors(10)=43.2743; // b5
    }
    else if (atom=="C")
    {
        descriptors( 0)= 6;     // Z
        descriptors( 1)= 0.0489; // a1
        descriptors( 2)= 0.2091; // a2
        descriptors( 3)= 0.7537; // a3
        descriptors( 4)= 1.1420; // a4
        descriptors( 5)= 0.3555; // a5
        descriptors( 6)= 0.1140; // b1
        descriptors( 7)= 1.0825; // b2
        descriptors( 8)= 5.4281; // b3
        descriptors( 9)=17.8811; // b4
        descriptors(10)=51.1341; // b5
    }
    else if (atom=="N")
    {
        descriptors( 0)= 7;     // Z
        descriptors( 1)= 0.0267; // a1
        descriptors( 2)= 0.1328; // a2
        descriptors( 3)= 0.5301; // a3
        descriptors( 4)= 1.1020; // a4
        descriptors( 5)= 0.4215; // a5
        descriptors( 6)= 0.0541; // b1
        descriptors( 7)= 0.5165; // b2
        descriptors( 8)= 2.8207; // b3
        descriptors( 9)=10.6297; // b4
        descriptors(10)=34.3764; // b5
    }
    else if (atom=="O")
    {
        descriptors( 0)= 8;     // Z
        descriptors( 1)= 0.0365; // a1
        descriptors( 2)= 0.1729; // a2
        descriptors( 3)= 0.5805; // a3
        descriptors( 4)= 0.8814; // a4
        descriptors( 5)= 0.3121; // a5
        descriptors( 6)= 0.0652; // b1
        descriptors( 7)= 0.6184; // b2
        descriptors( 8)= 2.9449; // b3
        descriptors( 9)= 9.6298; // b4
        descriptors(10)=28.2194; // b5
    }
    else if (atom=="P")
    {
        descriptors( 0)=15;     // Z
        descriptors( 1)= 0.1005; // a1
        descriptors( 2)= 0.4615; // a2
        descriptors( 3)= 1.0663; // a3
        descriptors( 4)= 2.5854; // a4
        descriptors( 5)= 1.2725; // a5
        descriptors( 6)= 0.0977; // b1
        descriptors( 7)= 0.9084; // b2
        descriptors( 8)= 4.9654; // b3
        descriptors( 9)=18.5471; // b4
        descriptors(10)=54.3648; // b5
    }
    else if (atom=="S")
    {
        descriptors( 0)=16;     // Z
        descriptors( 1)= 0.0915; // a1
        descriptors( 2)= 0.4312; // a2
        descriptors( 3)= 1.0847; // a3
        descriptors( 4)= 2.4671; // a4
        descriptors( 5)= 1.0852; // a5
        descriptors( 6)= 0.0838; // b1
        descriptors( 7)= 0.7788; // b2
        descriptors( 8)= 4.3462; // b3
        descriptors( 9)=15.5846; // b4
        descriptors(10)=44.6365; // b5
    }
    else if (atom=="Fe")
    {
        descriptors( 0)=26;     // Z
        descriptors( 1)= 0.1929; // a1
        descriptors( 2)= 0.8239; // a2
        descriptors( 3)= 1.8689; // a3
        descriptors( 4)= 2.3694; // a4
        descriptors( 5)= 1.9060; // a5
        descriptors( 6)= 0.1087; // b1
        descriptors( 7)= 1.0806; // b2
        descriptors( 8)= 4.7637; // b3
        descriptors( 9)=22.8500; // b4
        descriptors(10)=76.7309; // b5
    }
    else if (atom=="K")
    {
        descriptors( 0)=19;     // Z
        descriptors( 1)= 0.2149; // a1
        descriptors( 2)= 0.8703; // a2
        descriptors( 3)= 2.4999; // a3
        descriptors( 4)= 2.3591; // a4
        descriptors( 5)= 3.0318; // a5
        descriptors( 6)= 0.1660; // b1
        descriptors( 7)= 1.6906; // b2
        descriptors( 8)= 8.7447; // b3
        descriptors( 9)=46.7825; // b4
        descriptors(10)=165.6923; // b5
    }
    else if (atom=="F")
    {
        descriptors( 0)=9;     // Z
        descriptors( 1)= 0.0382; // a1
        descriptors( 2)= 0.1822; // a2
        descriptors( 3)= 0.5972; // a3
        descriptors( 4)= 0.7707; // a4
        descriptors( 5)= 0.2130; // a5
        descriptors( 6)= 0.0613; // b1
        descriptors( 7)= 0.5753; // b2
        descriptors( 8)= 2.6858; // b3
        descriptors( 9)= 8.8214; // b4
        descriptors(10)=25.6668; // b5
    }
    else if (atom=="Mg")
    {
        descriptors( 0)=12;     // Z
        descriptors( 1)= 0.1130; // a1
        descriptors( 2)= 0.5575; // a2
        descriptors( 3)= 0.9046; // a3
        descriptors( 4)= 2.1580; // a4
        descriptors( 5)= 1.4735; // a5
        descriptors( 6)= 0.1356; // b1
        descriptors( 7)= 1.3579; // b2
        descriptors( 8)= 6.9255; // b3
        descriptors( 9)=32.3165; // b4
        descriptors(10)=92.1138; // b5
    }
    else if (atom=="Cl")
    {
        descriptors( 0)=17;     // Z
        descriptors( 1)= 0.0799; // a1
        descriptors( 2)= 0.3891; // a2
        descriptors( 3)= 1.0037; // a3
        descriptors( 4)= 2.3332; // a4
        descriptors( 5)= 1.0507; // a5
        descriptors( 6)= 0.0694; // b1
        descriptors( 7)= 0.6443; // b2
        descriptors( 8)= 3.5351; // b3
        descriptors( 9)=12.5058; // b4
        descriptors(10)=35.8633; // b5
    }
    else if (atom=="Ca")
    {
        descriptors( 0)=20;     // Z
        descriptors( 1)= 0.2355; // a1
        descriptors( 2)= 0.9916; // a2
        descriptors( 3)= 2.3959; // a3
        descriptors( 4)= 3.7252; // a4
        descriptors( 5)= 2.5647; // a5
        descriptors( 6)= 0.1742; // b1
        descriptors( 7)= 1.8329; // b2
        descriptors( 8)= 8.8407; // b3
        descriptors( 9)=47.4583; // b4
        descriptors(10)=134.9613; // b5
    }
    else
        REPORT_ERROR(ERR_VALUE_INCORRECT,(std::string)"atomDescriptors: Unknown atom "+atom);
}

atomProjectionRadialProfile()

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

Atom projection radial profile. Returns the radial profile of the atom described by its profileCoefficients (Bspline coefficients).

atomRadialProfile()

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

Atom radial profile. Returns the radial profile of a given atom, i.e., the electron scattering factor convolved with a suitable low pass filter for sampling the volume at a sampling rate M*T. The radial profile is sampled at T Angstroms/pixel.

atomRadius()

double atomRadius

(

const std::string &

atom

)

Returns the radius of an atom. Returns 0 if the atom is not within the short list (H, C, N, O, S, P, Fe) of valid atoms.

The radius data is taken from http://www.webelements.com as the empirical radius.

Definition at line 168 of file pdb.cpp.

{
    switch (atom[0])
    {
    case 'H':
        return 0.25;
        break;
    case 'C':
        return 0.70;
        break;
    case 'N':
        return 0.65;
        break;
    case 'O':
        return 0.60;
        break;
    case 'P':
        return 1.00;
        break;
    case 'S':
        return 1.00;
        break;
    case 'E': // Iron Fe
        return 1.40;
        break;
    case 'K':
        return 2.20;
        break;
    case 'F':
        return 0.50;
        break;
    case 'G': // Magnesium Mg
        return 1.50;
        break;
    case 'L': // Chlorine Cl
        return 1.00;
        break;
    case 'A': // Calcium Ca
        return 1.80;
        break;
    default:
        return 0;
    }
}

computePDBgeometry()

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

Compute the center of mass and limits of a PDB file. The intensity column is used only for the pseudoatoms. It specifies from which column we should read the intensity. Valid columns are Bfactor or occupancy.

Definition at line 238 of file pdb.cpp.

{
    // Initialization
    centerOfMass.initZeros(3);
    limit0.resizeNoCopy(3);
    limitF.resizeNoCopy(3);
    limit0.initConstant(1e30);
    limitF.initConstant(-1e30);
    double total_mass = 0;

    // Initialize PDBRichPhantom and read atom struct file
    PDBRichPhantom pdbFile;

    // Read centered pdb
    pdbFile.read(fnPDB);

    // For each atom, correct necessary info
    bool useBFactor = intensityColumn=="Bfactor";
    for (auto& atom : pdbFile.atomList) {
        // Update center of mass and limits
        XX(limit0) = std::min(XX(limit0), atom.x);
        YY(limit0) = std::min(YY(limit0), atom.y);
        ZZ(limit0) = std::min(ZZ(limit0), atom.z);

        XX(limitF) = std::max(XX(limitF), atom.x);
        YY(limitF) = std::max(YY(limitF), atom.y);
        ZZ(limitF) = std::max(ZZ(limitF), atom.z);

        double weight;
        if (atom.name == "EN")
        {
            if (useBFactor)
                weight = atom.bfactor;
            else
                weight = atom.occupancy;
        }
        else
        {
            if (atom.record == "HETATM")
                continue;
            weight = (double) atomCharge(atom.name);
        }
        total_mass += weight;
        XX(centerOfMass) += weight * atom.x;
        YY(centerOfMass) += weight * atom.y;
        ZZ(centerOfMass) += weight * atom.z;
    }

    // Finish calculations
    centerOfMass /= total_mass;
}

distanceHistogramPDB()

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

Compute distance histogram of a PDB phantom. Consider the distance between each atom and its N nearest neighbours. Then, compute the histogram of these distances with Nbin samples.

Definition at line 1627 of file pdb.cpp.

{
    // Compute the histogram of distances
    const std::vector<Atom> &atoms=phantomPDB.atomList;
    int Natoms=atoms.size();
    MultidimArray<double> NnearestDistances;
    NnearestDistances.resize((Natoms-1)*Nnearest);
    for (int i=0; i<Natoms; i++)
    {
        std::vector<double> NnearestToThisAtom;
        const Atom& atom_i=atoms[i];
        for (int j=i+1; j<Natoms; j++)
        {
            const Atom& atom_j=atoms[j];
            double diffx=atom_i.x-atom_j.x;
            double diffy=atom_i.y-atom_j.y;
            double diffz=atom_i.z-atom_j.z;
            double dist=sqrt(diffx*diffx+diffy*diffy+diffz*diffz);
            if (maxDistance>0 && dist>maxDistance)
                continue;
            //std::cout << "Analyzing " << i << " and " << j << " -> d=" << dist << std::endl;
            size_t nearestSoFar=NnearestToThisAtom.size();
            if (nearestSoFar==0)
            {
                NnearestToThisAtom.push_back(dist);
                //std::cout << "Pushing d" << std::endl;
            }
            else
            {
                size_t idx=0;
                while (idx<nearestSoFar && NnearestToThisAtom[idx]<dist)
                    idx++;
                if (idx<nearestSoFar)
                {
                    NnearestToThisAtom.insert(NnearestToThisAtom.begin()+idx,1,dist);
                    if (NnearestToThisAtom.size()>Nnearest)
                        NnearestToThisAtom.erase(NnearestToThisAtom.begin()+Nnearest);
                }
                if (idx==nearestSoFar && nearestSoFar<Nnearest)
                {
                    NnearestToThisAtom.push_back(dist);
                    //std::cout << "Pushing d" << std::endl;
                }
            }
        }
        if (i<Natoms-1)
            for (size_t k=0; k<Nnearest; k++)
                NnearestDistances(i*Nnearest+k)=NnearestToThisAtom[k];
    }
    compute_hist(NnearestDistances, hist, 0, NnearestDistances.computeMax(), Nbins);
}

electronFormFactorFourier()

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

Compute the electron Form Factor in Fourier space. The electron scattering factor at a frequency f (Angstroms^-1) is computed as f_el(f)=sum_i(ai exp(-bi*x^2)).

Definition at line 1059 of file pdb.cpp.

{
    double retval=0;
    for (int i=1; i<=5; i++)
    {
        double ai=VEC_ELEM(descriptors,i);
        double bi=VEC_ELEM(descriptors,i+5);
        retval+=ai*exp(-bi*f*f);
    }
    return retval;
}

electronFormFactorRealSpace()

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

Compute the electron Form Factor in Real space. Konwing the electron form factor in Fourier space is easy to make an inverse Fourier transform and express it in real space. r is the distance to the center of the atom in Angstroms.

Definition at line 1087 of file pdb.cpp.

{
    double retval=0;
    for (int i=1; i<=5; i++)
    {
        double ai=descriptors(i);
        double bi=descriptors(i+5);
        double b=bi/(4*PI*PI);
        retval+=ai*sqrt(PI/b)*exp(-r*r/(4*b));
    }
    retval/=2*PI;
    return retval;
}

projectPDB()

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

Project PDB. Project the PDB following a certain projection direction.

Definition at line 1603 of file pdb.cpp.

{
    // Initialise projection
    proj().initZeros(Ydim, Xdim);
    proj().setXmippOrigin();
    proj.setAngles(rot, tilt, psi);

    // Compute volume to Projection matrix
    Matrix2D<double> VP = proj.euler;
    Matrix2D<double> PV = VP.inv();

    // Project all elements
    for (size_t i = 0; i < phantomPDB.getNumberOfAtoms(); i++)
    {
        try
        {
            projectAtom(phantomPDB.getAtom(i), proj, VP, PV, interpolator);
        }
        catch (XmippError &XE) {}
    }
}