src/common/sasa.cc

#include <stdio.h>
#include <math.h>
#include <string.h>
#include "open_prospect.h"
#include "sm_matrix.h"


float SasaCalc::CalcSasa( TargetStruct *target_data ) {
        float total = 0.0;
        if ( target_data->structure == NULL )
                return 0;
        
        for ( int i = 0; i < target_data->len; i++ ) {
                total += CalcResidueSasa( target_data, i );
        }
        return total;
}

PopsRSasa::PopsRSasa(ProspectParam *in_parent) : PopsSasa( in_parent, false ) {  
        char buffer[2000];
        params = parent->GetPopsRParam();
}


PopsSasa::PopsSasa(ProspectParam *in_parent, bool load)  : SasaCalc(), ParamElement( in_parent ) {  
        char buffer[2000];
        params = NULL;
        if (load) {
                params = parent->GetPopsParam();
        }
}



float PopsSasa::FindRadius( int atom_num ) {
        for ( int i = 0; i < params->param_count; i++) {
                if ( params->params[i].atom_num == atom_num ) {
                        return params->params[i].radius;
                }
        }
        return -1;
}

float PopsSasa::FindP( int atom_num ) {
        for ( int i = 0; i < params->param_count; i++) {
                if ( params->params[i].atom_num == atom_num ) {
                        return params->params[i].p_i;
                }
        }
        return -1;
}


#define PI 3.141592653589793
#define rSolv 1.6
float PopsSasa::CalcResidueSasa( TargetStruct *target_data, int i ) {
        float total = 0.0;      
        for ( int ia = 0; ia < target_data->structure[i].atom_count; ia++) {
                float ia_radius = FindRadius( target_data->structure[i].atom_name[ia] );
                if ( ia_radius > 0 ) {
                        float p_i = FindP(  target_data->structure[i].atom_name[ia] );                                  
                        float si = 4 * PI * pow( ia_radius + rSolv, 2 );                                        
                        float i_total = 1.0;
                        for ( int j = 0; j < target_data->len; j++ ) {
                                if ( i != j ) {
                                        for ( int ja = 0; ja < target_data->structure[j].atom_count; ja++) {
                                                float ja_radius = FindRadius( target_data->structure[j].atom_name[ja] );
                                                if ( ja_radius > 0 ) {
                                                        float r_ij = smDistance( target_data->structure[i].atom[ia], 
                                                                                                         target_data->structure[j].atom[ja] );
                                                        float b = 0;
                                                        if ( r_ij <= ia_radius + ja_radius + 2 * rSolv ) {
                                                                b = PI * ( ia_radius + rSolv ) * ( ia_radius + ja_radius + 2*rSolv - r_ij) *
                                                                (1 + ((ja_radius - ia_radius) / (r_ij)) );
                                                        }
                                                        i_total = i_total * (1 - ( p_i * params->p_14 * b) / si);
                                                }
                                        }               
                                        
                                }
                        }
                        total += si * i_total;
                }
        }
        return total;
}



float SimpleSasa::CalcResidueSasa( TargetStruct *target_data, int res ) {       
        return 0;       
}

float *ParseDsspAcc( FILE *file, int *count = NULL ) {
        char buffer[2000];
        char read_state = 0;
        int read_count = 0;     
        float *tmp_acc = new float[2000];
        while (fgets(buffer, sizeof(buffer), file)) {
                if ( !strncmp("  #", buffer, 3) ) {
                        read_state = 1;
                } else {
                        if ( read_state ) {
                                char tmp[ 5 ];
                                tmp[0] = buffer[35];
                                tmp[1] = buffer[36];
                                tmp[2] = buffer[37];
                                tmp[3] = 0;
                                tmp_acc[ read_count ] = atoi( tmp );
                                read_count++;
                        }                       
                }
        }
        float *acc_val = new float [ read_count ];
        for ( int i = 0; i < read_count; i++ ) {
                acc_val[i] = tmp_acc[i];                
        }
        delete tmp_acc;
        if ( count != NULL ) {
                (*count) =  read_count;
        }
        return acc_val;
}

void DsspSasa::LoadDssp(char *path) {           
        FILE *file = fopen( path, "r" );
        acc_val = ParseDsspAcc( file, &len );
        fclose( file );
}

float DsspSasa::CalcResidueSasa( TargetStruct *target_data, int res) { 
        if ( res >= len )
                return 0;
        return acc_val[ res ];
        
        
        /*
        char *text = target_data->WritePDBText( );
        
        char buffer[ 5000 ];
        
        sprintf(buffer, "dsspcmbi ");
        
        int pipe_to_dssp[2], dssp_from_dssp[2] ;
        pipe( pipe_to_dssp );
        pipe( pipe_from_dssp );

        if ( !(child_pid = fork()) ) {
                //this is the child, 
                
                //setup the pipe into dssp
                close(pipe_to_dssp[0]); //close the read end of the pipe
                close( stdin );
                dup2(pipe_to_dssp[1], stdin );
                
                char dssp_cmd[] = 
                dsspcmbi --
                
                
                write(pipe_files[1], xmltext, strlen(xmltext) );
                close(pipe_files[1]);
                exit(0);
        } else {
                close(pipe_files[1]); //close the write end of the pipe
                do {
                        read_file_and_append(pipe_files[0], &xmltext, &xmlsize);
                        status = wait4(child_pid, NULL, WNOHANG | WUNTRACED , NULL);
                        if (status == -1) {
                                fprintf(stderr, "Error while forking\n");
                                return NULL;
                        }
                        if (!status) 
                                usleep(300);
                        
                } while (status == 0 );
                read_file_and_append(pipe_files[0], &xmltext, &xmlsize);
                close(pipe_files[0]);
        }
        
        
        return acc_val[res];
         */
        
};

---