#include "pdb.h" double two_pi = M_PI + M_PI ; /* Insert an arc into the arcs array. each arc is two angles which are moved into the range -pi to pi, and arcsnum is the number of arcs. If the arc go over pi, it is cut to two arcs. Before actually inserting we check if it overlaps with any existing arc, and if it is extend the other arc to include both, remove it from arcs and reinsert it retunrn the number of arcs after the addition. */ insert_arc (double *arcs, double start, double end, int arcsnum) { if (start < (- M_PI - 0.00000001)) start += two_pi; if (end > M_PI + 0.00000001) end -= two_pi; if (start > end) { arcsnum = insert_arc(arcs, (- M_PI) , end,arcsnum) ; return insert_arc(arcs, start, M_PI, arcsnum) ; } else { int ii ; for(ii=0; ii< arcsnum; ii += 2) { double e_start = arcs[ii] ; double e_end = arcs[ii+1] ; if(start > e_start) { if (start <= e_end) { if (end > e_end) { arcs[ii] = arcs[ arcsnum-2]; arcs[ii+1] = arcs[ arcsnum-1]; return insert_arc(arcs,e_start,end,arcsnum-2) ; } else return arcsnum ; } } else if (end >= e_start) { if (end < e_end) end = e_end; arcs[ii] = arcs[ arcsnum-2]; arcs[ii+1] = arcs[ arcsnum-1]; return insert_arc(arcs,start,end,arcsnum-2) ; } } arcs[arcsnum++] = start; arcs[arcsnum++] = end; return arcsnum ; } } double sum_of_arcs(double *arcs, int arcsnum) { int ii = 0; double sum = 0.0 ; for (ii = 0 ; ii< arcsnum ; ii += 2) sum+= arcs[ii+1] - arcs[ii] ; return sum; } double zstep = 0.25 ; /* This setup all the relations betwen atom and its neibours in the xy plane, in parallel arrays. */ void setup_plane_data (atom_record *atom, atom_record **neibours, int neibour_count, double *xydists, double *xydists_sq, double *betas) { int ii; double xx = atom->x ; double yy = atom->y ; double zz = atom->z ; double deltax, deltay ,beta,sq; for (ii = 0 ; ii < neibour_count ; ii++) { atom_record *neibour = neibours [ ii] ; deltax = neibour->x - xx; deltay = neibour->y - yy; beta = atan2(deltay,deltax) ; betas[ii] = beta ; sq = deltax * deltax + deltay * deltay ; xydists[ii] = sqrt(sq) ; xydists_sq [ii] = sq; } } /* This is called on an atom, with a set of neibours and the count of them. It calculates the surface of an atom by: 1. Iterating over set of planes perpendicular to the Z axis, xy planes, seperated by zstep. For each plane: a. calculate the circle where the plane cut the sphere around the atom (plane_radius). the sphere is of extended radius, i.e. include the probe radius. for each neibour: calculate the cutting circle as above (neibour_plane_radius). if the circles are too far do nothing. if the central atom circle is around the other circle do nothing. if the central atom is inside the other circle, skip the plane. otherwise find the crssing points, by calculatin the angle between the line between the circles centres and one of the crossing points (alpha, by cosine law), and store the arc in arcs. sum the arcs in arcs, and substract from two pi to get the length of arc that is accessible for this atom in this plane, and add it to the total sum. 2. Multiply the total arc by zstep and the exetnded radius, to get the area accessible on the extended radius sphere, and then scale to the vdw radius. */ void calculate_atom_surface (atom_record *atom, atom_record **neibours, int neibour_count) { double xx = atom->x ; double yy = atom->y ; double zz = atom->z ; double xydists[200] ; double betas[200] ; double xydists_sq[200] ; atom_definition *ad= atom->definition ; double eradius = ad->eradius; double eradsq = ad->eradsq; double total_arc = 0.0; int ii,jj ; int grid_points_number = (int) (((eradius + eradius)/ zstep) + 0.5) ; double zgrid = zz - eradius - zstep *0.5 ; int arcsnum ; double arcs[200] ; setup_plane_data(atom,neibours,neibour_count, xydists, xydists_sq, betas) ; for (ii = 0 ; ii< grid_points_number ; ii++) { zgrid = zgrid + zstep ; { double aa =(zgrid - zz) ; double plane_radius_sq = eradsq - aa * aa ; double plane_radius ; arcsnum = 0 ; if (plane_radius_sq < 0.0) plane_radius_sq = 0.000001 ; plane_radius = sqrt(plane_radius_sq) ; for (jj = 0 ; jj< neibour_count ; jj++) { atom_record *neibour = neibours [jj] ; atom_definition * neibour_defn = neibour->definition ; double neibour_plane_radius ; double aa = (zgrid - neibour->z) ; double neibour_plane_radius_sq = neibour_defn->eradsq - aa * aa ; if (neibour_plane_radius_sq > 0.0) { double neibour_plane_radius = sqrt(neibour_plane_radius_sq) ; double xydist = xydists[ jj] ; if (xydist < neibour_plane_radius + plane_radius) { double temp = plane_radius - neibour_plane_radius ; if (xydist <= fabs(temp)) { if (temp > 0) continue; else goto next_plane; } else { double alpha = acos((xydists_sq[jj] + plane_radius_sq - neibour_plane_radius_sq) / (2.0 * xydist * plane_radius)) ; double beta = betas[jj] ; arcsnum= insert_arc(arcs,beta - alpha, beta+ alpha,arcsnum); } } } } total_arc += two_pi - sum_of_arcs(arcs,arcsnum) ; next_plane: ; } } atom->surface = total_arc * zstep * eradius ; } void call_calculate_atom_surface (atom_record *atom, double edges[3],int index) { if (! (atom->flags &( IGNORE_ATOM | DONT_CALCULATE_ATOM))) call_with_neibours(atom,calculate_atom_surface, 0) ; } extern double maximum_r ; void calculate_surface (file_records *file , double probe) { setup_file_records(file, probe,1) ; setup_minmax_xyzr(file,1) ; { double edge = maximum_r + maximum_r ; setup_neibouring_cube(file ,edge); traverse_file_atom_records (file, call_with_neibours, calculate_atom_surface) ; } }