/* EGAD: born_radius_sasa.cpp Navin Pokala and Tracy Handel Dept. of Molecular and Cell Biology University of California, Berkeley Copyright (C) 2003 Regents of the University of California GNU Public License Aug 12 2003 Absolutely no warranties are made or are implied with the use of this program or its parts. This file contains the born_radius_sasa function. */ #include "born_radius_sasa.h" /* calculates born radii and sasa for each atom in pdb (array of mini_pdbATOM) places these values in pdb[i].born_radius and pdb[i].sasa respectively */ void born_radius_sasa(mini_pdbATOM *pdb) { double r2; /* squared distance between two atoms */ double r; double xmax, ymax, zmax, xmin, ymin, zmin; int flag2, i_ctr, j_ctr; int i_heavy, j_heavy; mini_pdbATOM *firstAtom, *dihedAtom, *twobondAtom, *contactAtom; double dx, dy, dz; int i, j, free_vector_ptr_flag, sasa_flag; CARTESIAN_VECTOR *i_j_vector, *dummy_vector; mini_pdbATOM *micro_ptr[2]; SASA_SUM sasa_sum; extern CARTESIAN_VECTOR ***VECTOR_PTR; extern SASA_SUM SASA_SUM_TEMP; firstAtom = NULL; contactAtom = NULL; twobondAtom = NULL; dihedAtom = NULL; i_j_vector = NULL; dummy_vector = NULL; /* prevent crashes due to CARTESIAN_VECTOR matrix being too big */ if(MAX_ATOMS <= REAL_MAX_ATOMS) sasa_flag =1; else sasa_flag =2; free_vector_ptr_flag = 0; if(sasa_flag==1) { /* if this is NULL, then it means that this matrix should be created and destroyed internally */ if(VECTOR_PTR==NULL) /* for applications that require many calls to this function, this matrix should be allocated elsewhere */ { /* and kept for the duration of the need; this matrix is the largest object in terms of memory */ free_vector_ptr_flag=1; VECTOR_PTR = (CARTESIAN_VECTOR ***)calloc(MAX_ATOMS, sizeof(CARTESIAN_VECTOR **)); for(i=0;inext = NULL; micro_ptr[0] = &pdb[1]; micro_ptr[1] = NULL; } i_ctr = 1; j_ctr=1; while(pdb[i_ctr].seq_position!=ENDFLAG) /* initialize born_radius of charged atoms w/ born_bonded */ { if(pdb[i_ctr].atom_ptr->charge!=0) pdb[i_ctr].born_radius = pdb[i_ctr].atom_ptr->born_bonded; ++i_ctr; } i_ctr = 1; i_heavy=0; while(pdb[i_ctr].seq_position!=ENDFLAG) /* i loop */ { if(pdb[i_ctr].atom_ptr->atom_ptr->atom_class!=ENDFLAG) /* no U */ { if(sasa_flag==1) if(pdb[i_ctr].atom_ptr->atom_ptr->sasa_radius>0) { ++i_heavy; pdb[i_ctr].box = ((int)(pdb[i_ctr].coord.x/SASA_DISTANCE_CUTOFF)+500)*1000000 +((int)(pdb[i_ctr].coord.y/SASA_DISTANCE_CUTOFF)+500)*1000 +(int)(pdb[i_ctr].coord.z/SASA_DISTANCE_CUTOFF)+500 -(pdb[i_ctr].coord.x<0)*1000000 -(pdb[i_ctr].coord.y<0)*1000 -(pdb[i_ctr].coord.z<0); } xmax=pdb[i_ctr].coord.x+FORCEFIELD_DISTANCE_CUTOFF; xmin=pdb[i_ctr].coord.x-FORCEFIELD_DISTANCE_CUTOFF; ymax=pdb[i_ctr].coord.y+FORCEFIELD_DISTANCE_CUTOFF; ymin=pdb[i_ctr].coord.y-FORCEFIELD_DISTANCE_CUTOFF; zmax=pdb[i_ctr].coord.z+FORCEFIELD_DISTANCE_CUTOFF; zmin=pdb[i_ctr].coord.z-FORCEFIELD_DISTANCE_CUTOFF; j_ctr= i_ctr+1; j_heavy = i_heavy; while(pdb[j_ctr].seq_position!=ENDFLAG) /* j loop */ { if(pdb[j_ctr].atom_ptr->atom_ptr->atom_class!=ENDFLAG) { if(pdb[j_ctr].atom_ptr->atom_ptr->sasa_radius>0) ++j_heavy; /* only measure distances that are within a cube of side 2*FORCEFIELD_DISTANCE_CUTOFF centered at pdb[j_ctr].coord */ if(pdb[j_ctr].coord.x<=xmax) if(pdb[j_ctr].coord.x>=xmin) if(pdb[j_ctr].coord.y<=ymax) if(pdb[j_ctr].coord.y>=ymin) if(pdb[j_ctr].coord.z<=zmax) if(pdb[j_ctr].coord.z>=zmin) { dx = pdb[j_ctr].coord.x - pdb[i_ctr].coord.x; dy = pdb[j_ctr].coord.y - pdb[i_ctr].coord.y; dz = pdb[j_ctr].coord.z - pdb[i_ctr].coord.z; r2 = dx*dx + dy*dy + dz*dz; if(r2>= TINY && r2<=FORCEFIELD_DISTANCE_CUTOFF_SQRD) /* FORCEFIELD_DISTANCE_CUTOFF cutoff */ { flag2=0; r=sqrt(r2); if(sasa_flag==1) if(pdb[i_ctr].atom_ptr->atom_ptr->sasa_radius>0) if(pdb[j_ctr].atom_ptr->atom_ptr->sasa_radius>0) { if(r2<=SASA_RADII_SQRD_MATRIX[pdb[i_ctr].atom_ptr->atom_ptr->atom_class][pdb[j_ctr].atom_ptr->atom_ptr->atom_class]) { VECTOR_PTR[i_heavy][j_heavy] = i_j_vector; VECTOR_PTR[j_heavy][i_heavy] = i_j_vector; i_j_vector->i = i_heavy; i_j_vector->j = j_heavy; i_j_vector->r2 = r2; i_j_vector->r = r; dx/=r; dy/=r; dz/=r; i_j_vector->i_theta = acos(dz); i_j_vector->i_phi = atan2(dx,dy); dx = -dx; dy= -dy; dz= -dz; i_j_vector->j_theta = acos(dz); i_j_vector->j_phi = atan2(dx,dy); i_j_vector->next = (CARTESIAN_VECTOR *)malloc(sizeof(CARTESIAN_VECTOR)); i_j_vector = i_j_vector->next; i_j_vector->next = NULL; } } flag2 = get_atom_connectivity(r2, &pdb[i_ctr], &pdb[j_ctr], &firstAtom, &contactAtom, &twobondAtom, &dihedAtom); if(pdb[i_ctr].atom_ptr->charge!=0) /* i charged */ { if(flag2!=0) if(strcmp(pdb[j_ctr].atom_ptr->atom_ptr->atomtype, "HU")!=0) /* HU does not influence the born radius */ pdb[i_ctr].born_radius += born_nonbond_calc(pdb[i_ctr].atom_ptr, pdb[j_ctr].atom_ptr, r2); if(pdb[j_ctr].atom_ptr->charge!=0) /* i and j charged */ { if(flag2!=0) if(strcmp(pdb[i_ctr].atom_ptr->atom_ptr->atomtype, "HU")!=0) /* HU does not influence the born radius */ pdb[j_ctr].born_radius += born_nonbond_calc(pdb[j_ctr].atom_ptr, pdb[i_ctr].atom_ptr, r2); } } else if(pdb[j_ctr].atom_ptr->charge!=0) /* only j charged */ { if(flag2!=0) { if(strcmp(pdb[i_ctr].atom_ptr->atom_ptr->atomtype, "HU")!=0) /* HU does not influence the born radius */ pdb[j_ctr].born_radius += born_nonbond_calc(pdb[j_ctr].atom_ptr, pdb[i_ctr].atom_ptr, r2); } } } } } ++j_ctr; } } ++i_ctr; } i_ctr = 1; while(pdb[i_ctr].seq_position!=ENDFLAG) { if(pdb[i_ctr].atom_ptr->charge!=0) { if(pdb[i_ctr].born_radius >= 0) pdb[i_ctr].born_radius = -1.0; pdb[i_ctr].born_radius = -HALF_COULOMB_CONST/pdb[i_ctr].born_radius; } ++i_ctr; } if(sasa_flag==1) { sasa_sum = precalcSASAcalc(micro_ptr); i_j_vector = dummy_vector; freeCARTESIAN_VECTOR(i_j_vector); for(i=1;i<=i_heavy;++i) for(j=1;j<=i_heavy;++j) { if(VECTOR_PTR[i][j]!=NULL) { VECTOR_PTR[i][j]=NULL; VECTOR_PTR[j][i]=NULL; } } if(free_vector_ptr_flag==1) { for(i=0;i