voronoi.h File Reference

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <sys/types.h>
#include <unistd.h>
#include "rpdb.h"
#include "writepdb.h"
#include "calc.h"
#include "utils.h"
#include "../src/qhull/qvoronoi.h"
#include "memhandler.h"

Go to the source code of this file.

Data Structures
struct	s_vvertice
struct	s_lst_vvertice
Defines
#define	M_VORONOI_SUCCESS   0
#define	M_APOLAR_AS   0
#define	M_POLAR_AS   1
#define	M_PREC_TOLERANCE   1e-5
#define	M_BUFSIZE   1e7
Functions
s_lst_vvertice *	load_vvertices (s_pdb *pdb, int min_apol_neigh, float ashape_min_size, float ashape_max_size)
float	testVvertice (float xyz[3], int curNbIdx[4], s_atm *atoms, float min_asph_size, float max_asph_size, s_lst_vvertice *lvvert)
void	set_barycenter (s_vvertice *v)
int	is_in_lst_vert (s_vvertice **lst_vert, int nb_vert, int v_id)
int	is_in_lst_vert_p (s_vvertice **lst_vert, int nb_vert, s_vvertice *vert)
void	write_pqr_vert (FILE *f, s_vvertice *v)
void	write_pdb_vert (FILE *f, s_vvertice *v)
float	get_verts_volume_ptr (s_vvertice **verts, int nvert, int niter, float correct)
void	print_vvertices (FILE *f, s_lst_vvertice *lvvert)
void	free_vert_lst (s_lst_vvertice *lvvert)

---

Define Documentation

#define M_APOLAR_AS   0

alpha sphere type - hydrophilic alpha sphere

Definition at line 57 of file voronoi.h.

Referenced by fill_vvertices(), get_vert_apolar_density(), set_descriptors(), updateIds(), write_pdb_vert(), and write_pqr_vert().

#define M_BUFSIZE   1e7

buffer size

Definition at line 61 of file voronoi.h.

#define M_POLAR_AS   1

tolerance for coordinate imprecion during alpha sphere search

Definition at line 58 of file voronoi.h.

Referenced by fill_vvertices().

#define M_PREC_TOLERANCE   1e-5

Definition at line 59 of file voronoi.h.

Referenced by testVvertice().

#define M_VORONOI_SUCCESS   0

alpha sphere type - hydrophobic alpha sphere

Definition at line 56 of file voronoi.h.

Referenced by load_vvertices().

---

Function Documentation

void free_vert_lst

(

s_lst_vvertice *

lvvert

)

## FUNCTION: free_vert_lst

## SPECIFICATION: Free memory

## PARAMETERS: @ s_lst_vvertice *lvvert : Data to free

## RETURN: void

Definition at line 530 of file voronoi.c.

References s_lst_vvertice::h_tr, my_free(), s_lst_vvertice::pvertices, s_lst_vvertice::tr, and s_lst_vvertice::vertices.

Referenced by c_lst_pocket_free().

{
    if(lvvert) {

        if(lvvert->vertices) {
            my_free(lvvert->vertices) ;
            lvvert->vertices = NULL ;
        }
        if(lvvert->pvertices) {
            my_free(lvvert->pvertices) ;
            lvvert->pvertices = NULL ;
        }
        if(lvvert->tr) {
                my_free(lvvert->tr) ;
                lvvert->tr = NULL ;
        }
        if(lvvert->h_tr) {
                my_free(lvvert->h_tr) ;
                lvvert->h_tr = NULL ;
        }
        my_free(lvvert) ;
    }
}

float get_verts_volume_ptr	(	s_vvertice **	verts,
		int	nvert,
		int	niter,
		float	correct
	)

## FUNCTION: get_verts_volume_ptr

## SPECIFICATION: Get an monte carlo approximation of the volume occupied by the alpha spheres given in argument (list of pointers)

## PARAMETRES: @ s_vvertice **verts: List of pointer to alpha spheres @ int nvert: Number of spheres @ int niter: Number of monte carlo iteration to perform @ float correct: radius for which the size of an alpha sphere should be corrected in order to calculate the volume

## RETURN: float: volume.

Definition at line 454 of file voronoi.c.

References rand_uniform(), s_vvertice::ray, s_vvertice::x, s_vvertice::y, and s_vvertice::z.

Referenced by set_descriptors().

{
        int i = 0, j = 0,
                nb_in = 0;

        float xmin = 0.0, xmax = 0.0,
                  ymin = 0.0, ymax = 0.0,
                  zmin = 0.0, zmax = 0.0,
                  xtmp = 0.0, ytmp = 0.0, ztmp = 0.0,
                  xr   = 0.0, yr   = 0.0, zr   = 0.0,
                  vbox = 0.0 ;

        s_vvertice *vcur = NULL ;

        /* First, search extrems coordinates to get a contour box of the molecule */
        for(i = 0 ; i < nvert ; i++) {
                vcur = verts[i] ;

                if(i == 0) {
                        xmin = vcur->x - vcur->ray + correct ; xmax = vcur->x + vcur->ray + correct ;
                        ymin = vcur->y - vcur->ray + correct ; ymax = vcur->y + vcur->ray + correct ;
                        zmin = vcur->z - vcur->ray + correct ; zmax = vcur->z + vcur->ray + correct ;
                }
                else {
                /* Update the minimum and maximum extreme point */
                        if(xmin > (xtmp = vcur->x - vcur->ray + correct)) xmin = xtmp ;
                        else if(xmax < (xtmp = vcur->x + vcur->ray + correct)) xmax = xtmp ;

                        if(ymin > (ytmp = vcur->y - vcur->ray + correct)) ymin = ytmp ;
                        else if(ymax < (ytmp = vcur->y + vcur->ray + correct)) ymax = ytmp ;

                        if(zmin > (ztmp = vcur->z - vcur->ray + correct)) zmin = ztmp ;
                        else if(zmax < (ztmp = vcur->z + vcur->ray + correct)) zmax = ztmp ;
                }
        }

        /* Next calculate the contour box volume */
        vbox = (xmax - xmin)*(ymax - ymin)*(zmax - zmin) ;

        /* Then apply monte carlo approximation of the volume.   */
        for(i = 0 ; i < niter ; i++) {
                xr = rand_uniform(xmin, xmax) ;
                yr = rand_uniform(ymin, ymax) ;
                zr = rand_uniform(zmin, zmax) ;

                for(j = 0 ; j < nvert ; j++) {
                        vcur = verts[j] ;
                        xtmp = vcur->x - xr ;
                        ytmp = vcur->y - yr ;
                        ztmp = vcur->z - zr ;

                /* Compare r^2 and dist(center, random_point)^2 */
                        if(((correct+vcur->ray)*(vcur->ray + correct)) > (xtmp*xtmp + ytmp*ytmp + ztmp*ztmp)) {
                        /* The point is inside one of the vertice!! */
                                nb_in ++ ; break ;
                        }
                }
        }

        /* Ok lets just return the volume Vpok = Nb_in/Niter*Vbox */
        return ((float)nb_in)/((float)niter)*vbox;
}

int is_in_lst_vert	(	s_vvertice **	lst_vert,
		int	nb_vert,
		int	v_id
	)

## FUNCTION: is_in_lst_vert

## SPECIFICATION: Says wether a vertice of id v_id is in a list of vertices or not

## PARAMETRES:

## RETURN: 1 if the vertice is in the tab, 0 if not

Definition at line 567 of file voronoi.c.

Referenced by get_mol_vert_neigh().

{
        int i ;
        for(i = 0 ; i < nb_vert ; i++) {
                if(v_id == lst_vert[i]->id) return 1 ;
        }

        return 0 ;
}

int is_in_lst_vert_p	(	s_vvertice **	lst_vert,
		int	nb_vert,
		s_vvertice *	vert
	)

## FUNCTION: is_in_lst_vert

## SPECIFICATION: Says wether a vertice of id v_id is in a list of vertices or not

## PARAMETRES:

## RETURN: 1 if the vertice is in the tab, 0 if not

Definition at line 590 of file voronoi.c.

{
        int i ;
        for(i = 0 ; i < nb_vert ; i++) {
                if(vert == lst_vert[i]) return 1 ;
        }

        return 0 ;
}

s_lst_vvertice* load_vvertices	(	s_pdb *	pdb,
		int	min_apol_neigh,
		float	asph_min_size,
		float	asph_max_size
	)

## FUNCTION: s_lst_vvertice

## SPECIFICATION: Calculate voronoi vertices using an ensemble of atoms, and then load resulting vertices into a s_lst_vvertice structure. The function call an external programm qvoronoi, part of qhull programme which can be download at: http://www.qhull.org/download/ or installed with apt-get install qhull-bin

## PARAMETRES: @ s_pdb *pdb : PDB informations @ int min_apol_neigh : Number of apolar neighbor of a vertice to be considered as apolar @ float asph_min_size : Minimum size of voronoi vertices to retain @ float asph_max_size : Maximum size of voronoi vertices to retain

## RETURN: s_lst_vvertice * :The structure containing the list of vertices.

Definition at line 90 of file voronoi.c.

References fill_vvertices(), s_lst_vvertice::h_tr, s_pdb::latoms, M_VORONOI_SUCCESS, my_free(), my_malloc(), my_realloc(), s_lst_vvertice::n_h_tr, s_pdb::natoms, s_atm::symbol, s_atm::x, s_atm::y, and s_atm::z.

Referenced by search_pocket().

{
        int i, nb_h=0;
        s_atm *ca = NULL ;
        s_lst_vvertice *lvvert = NULL ;
        char tmpn1[250]="";
        char tmpn2[250]="";
        pid_t pid=getpid();
        sprintf(tmpn1,"/tmp/qvoro_tmp_%d.dat",pid);
        sprintf(tmpn2,"/tmp/fpocket_qvor_%d.dat",pid);
        FILE *fvoro = fopen(tmpn1, "w+");
        FILE *ftmp=fopen(tmpn2,"w");
        /*FILE *ftmp=fopen("/tmp/fpocket_qvor.dat","w");
        FILE *fvoro = fopen("/tmp/qvoro_tmp.dat", "w+");*/
/*      lvvert->h_tr=(int *)my_malloc(sizeof(int));*/



        if(fvoro != NULL) {
                lvvert = (s_lst_vvertice *)my_malloc(sizeof(s_lst_vvertice)) ;
                lvvert->h_tr=NULL;
                /* Loop a first time to get out how many heavy atoms are in the file */
                for(i = 0; i <  pdb->natoms ; i++){
                        ca = (pdb->latoms)+i ;
                        if(strcmp(ca->symbol,"H")) {
                                lvvert->h_tr=(int *)my_realloc(lvvert->h_tr,sizeof(int)*(i-nb_h+1)) ;
                                lvvert->h_tr[i-nb_h]=i ;
                        }
                        else nb_h++;
                }
                lvvert->n_h_tr=i-nb_h;

                /* Write the header for qvoronoi */
                fprintf(fvoro,"3 rbox D3\n%d\n", lvvert->n_h_tr);
                /* Loop a second time for the qvoronoi input coordinates */
                for(i = 0; i <  pdb->natoms ; i++){
                        ca = (pdb->latoms)+i ;
                        if(strcmp(ca->symbol,"H")) {
                        /* Only if this is a heavy atom export it for voronoi tesselation,
                         * else discard it */
                                fprintf(fvoro,"%.3f %.3f %.3f \n", ca->x, ca->y, ca->z);
                        }
                }

                fflush(fvoro) ;
                rewind(fvoro);
                
                //int status = system("qvoronoi p i Pp Fn < voro_tmp.dat > voro.tmp") ;
                run_qvoronoi(fvoro,ftmp);
                int status=M_VORONOI_SUCCESS;


              
                if(status == M_VORONOI_SUCCESS) {
                        fill_vvertices(lvvert, tmpn2, pdb->latoms, pdb->natoms,
                                                        min_apol_neigh, asph_min_size, asph_max_size);
                }
                else {
                        my_free(lvvert);
                        lvvert = NULL ;
                        fprintf(stderr, "! Voronoi command failed with status %d...\n", status) ;
                }
                
          
        }
        else {
                fprintf(stderr, "! File for Voronoi vertices calculation couldn't be opened...\n") ;
        }
        fclose(fvoro) ;
        fclose(ftmp);
        remove(tmpn1);
        remove(tmpn2);
        return lvvert ;
}

void print_vvertices	(	FILE *	f,
		s_lst_vvertice *	lvvert
	)

## FUNCTION: print_vvertices

## SPECIFICATION: Print function.

## PARAMETERS: @ FILE *f : Buffer to print in @ s_lst_vvertice *lvvert : Vertices to print

## RETURN: void

Definition at line 411 of file voronoi.c.

References dist(), s_atm::id, s_vvertice::neigh, s_lst_vvertice::nvert, s_vvertice::sort_x, s_lst_vvertice::vertices, s_atm::x, s_vvertice::x, s_atm::y, s_vvertice::y, s_atm::z, and s_vvertice::z.

{
        if(lvvert) {
                if(lvvert->vertices) {
                        int i ;
                        for(i = 0 ; i < lvvert->nvert ; i++) {
                                s_vvertice *v = &(lvvert->vertices[i]) ;
                                if( v->neigh[0] &&  v->neigh[1] && v->neigh[2] &&  v->neigh[3]) {
                                        fprintf(f, "====== Vertice %d: =====\n", i);
                                        fprintf(f, "- x = %f, y = %f, z = %f\n", v->x, v->y, v->z);
                                        fprintf(f, "- ix = %d\n",v->sort_x);

                                        float d1 = dist(v->x, v->y, v->z, v->neigh[0]->x, v->neigh[0]->y, v->neigh[0]->z) ;
                                        float d2 = dist(v->x, v->y, v->z, v->neigh[1]->x, v->neigh[1]->y, v->neigh[1]->z) ;
                                        float d3 = dist(v->x, v->y, v->z, v->neigh[2]->x, v->neigh[2]->y, v->neigh[2]->z) ;
                                        float d4 = dist(v->x, v->y, v->z, v->neigh[3]->x, v->neigh[3]->y, v->neigh[3]->z) ;

                                        fprintf(f, "- Neighbour: \n1 - %f (%f %f %f: %d) \n2 - %f (%f %f %f: %d)\n3 - %f (%f %f %f: %d)\n4 - %f (%f %f %f: %d)\n", d1, v->neigh[0]->x, v->neigh[0]->y, v->neigh[0]->z, v->neigh[0]->id, d2, v->neigh[1]->x, v->neigh[1]->y, v->neigh[1]->z, v->neigh[1]->id, d3, v->neigh[2]->x,  v->neigh[2]->y ,  v->neigh[2]->z, v->neigh[2]->id, d4, v->neigh[3]->x, v->neigh[3]->y, v->neigh[3]->z, v->neigh[3]->id);
                                }
                        }
                }
        }
}

void set_barycenter

(

s_vvertice *

v

)

## FUNCTION: set_barycenter

## SPECIFICATION: Set barycenter of a vertice using it's 4 contacting atoms.

## PARAMETERS: @ s_vvertice *v: The vertice

## RETURN: void

Definition at line 321 of file voronoi.c.

References s_vvertice::bary, s_vvertice::neigh, s_atm::x, s_atm::y, and s_atm::z.

Referenced by fill_vvertices().

{
        int i ;
        float xsum = 0.0,
                  ysum = 0.0,
                  zsum = 0.0 ;

        for(i = 0 ; i < 4 ; i++) {
                xsum += v->neigh[i]->x ;
                ysum += v->neigh[i]->y ;
                zsum += v->neigh[i]->z ;
        }

        v->bary[0] = xsum*0.25 ;
        v->bary[1] = ysum*0.25 ;
        v->bary[2] = zsum*0.25 ;

}

float testVvertice	(	float	xyz[3],
		int	curNbIdx[4],
		s_atm *	atoms,
		float	min_asph_size,
		float	max_asph_size,
		s_lst_vvertice *	lvvert
	)

## FUNCTION: testVvertice

## SPECIFICATION: Test if alpha sphere conditions are fulfilled for current vertice

## PARAMETERS: @ float xyz[3] : Coordinates of current vertice @ int curNbIdx[4] : Indexes of atomic neighbours of the current vertice @ s_atm *atoms : List of all atoms @ float min_asph_size : Minimum size of alpha spheres. @ float max_asph_size : Maximum size of alpha spheres.

## RETURN: float : -1 if conditions are not fulfilled, else the alpha sphere radius is returned.

Definition at line 359 of file voronoi.c.

References dist(), s_lst_vvertice::h_tr, M_PREC_TOLERANCE, s_atm::x, s_atm::y, and s_atm::z.

Referenced by fill_vvertices().

{
        float x = xyz[0],
                  y = xyz[1],
                  z = xyz[2] ;

        s_atm *cura = &(atoms[lvvert->h_tr[curNbIdx[0]]]) ;

        float distVatom1 = dist(x, y, z, cura->x, cura->y, cura->z) ;
        float distVatom2,
                  distVatom3,
                  distVatom4;

        if(min_asph_size <= distVatom1 + M_PREC_TOLERANCE  && distVatom1 - M_PREC_TOLERANCE <= max_asph_size){
                cura = &(atoms[lvvert->h_tr[curNbIdx[1]]]) ;
                distVatom2 = dist(x, y, z, cura->x, cura->y, cura->z);

                cura = &(atoms[lvvert->h_tr[curNbIdx[2]]]) ;
                distVatom3 = dist(x, y, z, cura->x, cura->y, cura->z);

                cura = &(atoms[lvvert->h_tr[curNbIdx[3]]]) ;
                distVatom4=dist(x, y, z, cura->x, cura->y, cura->z);

                /* Test if all 4 neighbours are on the alpha sphere surface
                 * (approximate test) */
                if(fabs(distVatom1-distVatom2) < M_PREC_TOLERANCE &&
                   fabs(distVatom1-distVatom3)  < M_PREC_TOLERANCE &&
                   fabs(distVatom1-distVatom4) < M_PREC_TOLERANCE){

                        return distVatom1;
                }

        }
        return -1.0;
}

void write_pdb_vert	(	FILE *	f,
		s_vvertice *	v
	)

-- --

OUTPUT FUNCTIONS

-- -- -- ## FUNCTION: void write_pdb_vertice(FILE *f, s_vvertice *v)

## SPECIFICATION: Write a voronoi vertice in pdb format.

## PARAMETRES: @ FILE *f: file to write in @ s_vvertice *v: The vertice

## RETURN:

Definition at line 624 of file voronoi.c.

References s_vvertice::id, M_APOLAR_AS, s_vvertice::resid, s_vvertice::type, write_pdb_atom_line(), s_vvertice::x, s_vvertice::y, and s_vvertice::z.

Referenced by write_pockets_single_pdb().

{
        if(v->type==M_APOLAR_AS) write_pdb_atom_line(f, "HETATM", v->id, "APOL",
                                                                                                 ' ', "STP", "C", v->resid, ' ',
                                                                                                 v->x, v->y, v->z, 0.0, 0.0,
                                                                                                 "Ve", -1);

        else write_pdb_atom_line(f, "HETATM", v->id, " POL", ' ', "STP", "C",
                                                                v->resid, ' ',v->x, v->y, v->z,0.0, 0.0,
                                                                "Ve", -1);
}

void write_pqr_vert	(	FILE *	f,
		s_vvertice *	v
	)

## FUNCTION: write_pqr_vertice

## SPECIFICATION: Write a voronoi vertice in pqr format.

## PARAMETRES: @ FILE *f : file to write in @ s_vvertice *v : The vertice

## RETURN: void

Definition at line 651 of file voronoi.c.

References s_vvertice::id, M_APOLAR_AS, s_vvertice::ray, s_vvertice::resid, s_vvertice::type, write_pqr_atom_line(), s_vvertice::x, s_vvertice::y, and s_vvertice::z.

Referenced by mdpocket_characterize(), write_mdpockets_concat_pqr(), write_pocket_pqr(), write_pocket_pqr_DB(), and write_pockets_single_pqr().

{
        if(v->type==M_APOLAR_AS) write_pqr_atom_line(f, "ATOM", v->id, "APOL", ' ',
                                                                                                 "STP", " ", v->resid, ' ',
                                                                                                  v->x, v->y, v->z, 0.0, v->ray);

        else write_pqr_atom_line(f, "ATOM", v->id, " POL", ' ', "STP", " ",
                                                         v->resid, ' ',v->x, v->y, v->z,0.0, v->ray);
}