#include "pdb.h" /* code to calculate and use a neibouring cube, that makes finding neibouring atons faster */ double cube_edges [3] ; void setup_atom_definition (atom_definition * ad,double *probe) { double eradius = ad->vdw + *probe ; ad->eradius = eradius ; ad->eradsq = eradius * eradius ; if (eradius > maximum_r) maximum_r = eradius ; if (eradius < minimum_r) minimum_r = eradius ; } void check_min_max_xyz (atom_record *atom, void *arg,int index) { if (! (atom->flags & IGNORE_ATOM)) { double y = atom->y ; double z = atom->z ; double x = atom->x ; if (x > maximum_x) maximum_x = x ; if (x < minimum_x) minimum_x = x ; if (y > maximum_y) maximum_y = y ; if (y < minimum_y) minimum_y = y ; if (z > maximum_z) maximum_z = z ; if (z < minimum_z) minimum_z = z ; } } /* this setup the minimum and maximum x/y/z/r, */ void setup_minmax_xyzr (file_records *file,int new ) { if (new) { minimum_x = minimum_y = minimum_z = 999999.0 ; maximum_x = maximum_y = maximum_z = -999999.0 ; } traverse_file_atom_records(file, check_min_max_xyz , 0) ; } /* and also the extended radius and its square in the atom-definitions. This needs to be called before using the neibouring cube. */ void setup_file_records (file_records *file, double probe,int new) { if(new) { minimum_r = 999999.0 ; maximum_r = -999999.0 ; } traverse_atom_definitions(setup_atom_definition , &probe) ; } #define CUBE_SIDE 30 atom_record *n_cube[CUBE_SIDE][CUBE_SIDE][CUBE_SIDE] ; /* this is called for each atom, and put it in n_cube in the right place. The atom chaining field is set to point to the exisiting entry, so in the end each each entry is pointing to chain of atoms */ void put_atom_in_cube(atom_record *atom, double edges[3]) { if (! (atom->flags & IGNORE_ATOM)) { int x =(int)floor( (atom->x - minimum_x) /edges[0]) ; int y =(int)floor( (atom->y - minimum_y) /edges[1]) ; int z =(int)floor( (atom->z - minimum_z) /edges[2]) ; atom->chaining = n_cube[x][y][z] ; if(debugging & DEBUG_CUBE) { print_an_atom(atom,stderr); printf("storing %d %d %d\n", x , y , z) ; } n_cube[x][y][z] = atom ; } } void add_file_to_neibouring_cube(file_records *file) { traverse_file_atom_records(file , put_atom_in_cube, cube_edges) ; } /* Setup a grid covering the muleculle, such that neibouring atoms of an atom in one of the small cubes are all in the neibouring cubes. this make finding all the neibours of an atom cheaper. The atoms in each small cube are chained using the chaining field in the atoms. the values of miimum/maximum of x/y/z have to already setup when this called, by setup_file_records above or similar call. the file_records contain the atoms, the egdes the desired edges for x y and z. If this are too small for the cube to include all the volume, the egdes are increased. the edges notmally should be used later by call_with_neibours */ void in_setup_neibouring_cube (file_records *file, double edge) { double deltax = maximum_x - minimum_x + 0.1; double deltay = maximum_y - minimum_y + 0.1; double deltaz = maximum_z - minimum_z + 0.1; double edgex = cube_edges[0] = edge ; double edgey = cube_edges[1] = edge ; double edgez = cube_edges[2] = edge ; { int ii,jj,kk ; for(ii= 0 ; ii< CUBE_SIDE ; ii++) for(jj= 0 ; jj< CUBE_SIDE ; jj++) for(kk= 0 ; kk< CUBE_SIDE ; kk++) n_cube[ii][jj][kk] = NULL ; } if (edgex < deltax / CUBE_SIDE) { edgex = cube_edges[0] = deltax/CUBE_SIDE ; if (debugging & 1 ) fprintf(stderr, LF_STRING_03, deltax, CUBE_SIDE,edgex); } if (edgey < deltay / CUBE_SIDE) { edgey = cube_edges[1] = deltay/CUBE_SIDE ; if (debugging & 1 ) fprintf(stderr, LF_STRING_04, deltay, CUBE_SIDE,edgey); } if (edgez < deltaz / CUBE_SIDE) { edgez = cube_edges[2] = deltaz/CUBE_SIDE ; if (debugging & 1 ) fprintf(stderr, LF_STRING_05, deltaz, CUBE_SIDE,edgez); } add_file_to_neibouring_cube(file) ; } void setup_neibouring_cube (file_records *file, double edge) { setup_minmax_xyzr(file,1) ; in_setup_neibouring_cube(file, edge) ; } double distance_check = -1.0 ; /* call the functions with the atom and an array containing all its neiboring. This is using the cube setup by setup_neibouring_cube */ void call_with_neibours(atom_record *atom,void (*function)() , int what ) { int ii = SELECTED_ATOM_P (atom) ; switch(what) { case CWN_ALL_ATOMS :ii = 1; break ; case CWN_NOT_IGNORE : ii = !(atom->flags & IGNORE_ATOM) ; break; case CWN_NOT_SELCTED : ii = ii ? 0 : 1 ; break ; } if (ii) { atom_record *neibours[1000] ; int neibour_count = 0 ; int x =(int)floor( (atom->x - minimum_x) /cube_edges[0]) ; int y =(int)floor( (atom->y - minimum_y) /cube_edges[1]) ; int z =(int)floor( (atom->z - minimum_z) /cube_edges[2]) ; int xmin = x -1 ; int xmax = x +1 ; int ymin = y -1 ; int ymax = y +1 ; int zmin = z -1 ; int zmax = z +1 ; int yy, zz; if (0 > xmin) xmin = 0 ; if (xmax >= CUBE_SIDE ) { xmax-- ; if (xmax >= CUBE_SIDE) xmax--; } if (0 > ymin) ymin = 0 ; if (ymax >= CUBE_SIDE ) { ymax-- ; if (ymax >= CUBE_SIDE ) ymax-- ; } if (0 > zmin) zmin = 0 ; if (zmax >= CUBE_SIDE ) { zmax-- ; if (zmax >= CUBE_SIDE ) zmax-- ; } for ( ; xmin <= xmax; xmin++) for (yy = ymin ; yy <= ymax; yy++) for ( zz= zmin ; zz <= zmax; zz++) { atom_record *neibour = n_cube[xmin][yy][zz] ; while(neibour != NULL) { if ((int)neibour != (int)atom && ((distance_check > 0.0 ? distance_check : (atom->definition->eradius + neibour->definition->eradius)) > atoms_distance(atom,neibour))) neibours[ neibour_count ++] = neibour ; neibour = neibour->chaining ; } } (*function)(atom,neibours, neibour_count) ; } } void always_call_with_neibours(atom_record *atom,void (*function)() , int what ) { call_with_neibours(atom, function, CWN_ALL_ATOMS) ; }