/* EGAD: energy_profile_table.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 functions for calculating and outputing energy profile tables. diagonal is internal and solvation energy for the residue of interest bkbn-X = interaction energy between bkbn and X res(X)-res(Y) = interaction energy between sidechains X and Y 0 1 2 n-1 n n+1 bkbn res(1) res(2) ..... res(n-1) res(n) totals 0 bkbn bkbn_total 1 res(1) res(1)_total 2 res(2) res(2)_total . . . n-1 res(n-1) res(n-1)_total n res(n) res(n)_total n+1 totals bkbn res(1) res(2) ..... res(n-1) res(n) total energy for the structure total total total total total */ #include "energy_profile_table.h" void energy_profile_table_difference(int num_res, double **energy_table_diff, double **energy_table_A, double **energy_table_B) { int i,j; if(energy_table_diff==NULL) failure_report("ERROR energy_table_diff must be allocated for energy_profile_table.cpp: energy_profile_table_difference", "abort"); for(i=0;i<=num_res+1;++i) for(j=0;j<=num_res+1;++j) energy_table_diff[i][j] = energy_table_A[i][j] - energy_table_B[i][j]; } /* output energy_table to file name.energy_profile_table; called by output_PROTEIN */ void output_energy_profile_table(double **energy_table, char *name, int num_res, SEQPOS_TEXT_MAPPING_LIST *seqpos_text_map, char *sequence) { FILE *output; char *filename; int i,j; if(energy_table==NULL) { fprintf(stderr,"ERROR energy_profile_table not created; cannot output_energy_profile_table\n"); fprintf(stderr,"Dumping core....\n"); abort(); } filename = (char *)calloc(MXLINE_INPUT,sizeof(char)); sprintf(filename, "%s.energy_profile_table", name); output = fopen_file(filename,"w"); /* labels */ fprintf(output,"\tbkbn\t"); for(i=1;i<=num_res;++i) fprintf(output,"%s %c\t",seqpos_text_map[i].seqpos_text, sequence[i-1]); fprintf(output,"total\n"); fprintf(output,"bkbn\t%f\t",energy_table[0][0]); for(i=1;i<=num_res;++i) fprintf(output,"%f\t",energy_table[0][i]); fprintf(output,"%f\n",energy_table[0][num_res+1]); for(i=1;i<=num_res;++i) { fprintf(output,"%s %c\t",seqpos_text_map[i].seqpos_text, sequence[i-1]); for(j=0;j<=num_res;++j) fprintf(output,"%f\t",energy_table[i][j]); fprintf(output,"%f\n",energy_table[i][num_res+1]); } fprintf(output,"total\t"); for(i=0;i<=num_res;++i) fprintf(output,"%f\t",energy_table[num_res+1][i]); fprintf(output,"%f\n",energy_table[num_res+1][num_res+1]); fclose(output); free_memory(filename); } /* place backbone atoms (not CB) from pdb into backbone array (allocated by calling function). for each sidechain i, place sidechain atoms (including CB) into sidechain[i] array. (**sidechain allocated by calling function). */ void extract_sub_structures(mini_pdbATOM *pdb, mini_pdbATOM *backbone, mini_pdbATOM **sidechain) { int i,b,c,index,k; i=1; b=1; index=0; while(pdb[i].seq_position!=ENDFLAG) { /* backbone atoms (not CB) */ if(strcmp(pdb[i].atom_ptr->atomname,"N")==0 || strcmp(pdb[i].atom_ptr->atomname,"H")==0 || strcmp(pdb[i].atom_ptr->atomname,"1H")==0 || strcmp(pdb[i].atom_ptr->atomname,"2H")==0 || strcmp(pdb[i].atom_ptr->atomname,"3H")==0 || strcmp(pdb[i].atom_ptr->atomname,"CA")==0 || strcmp(pdb[i].atom_ptr->atomname,"HA")==0 || strcmp(pdb[i].atom_ptr->atomname,"C")==0 || strcmp(pdb[i].atom_ptr->atomname,"O")==0 || strcmp(pdb[i].atom_ptr->atomname,"OT")==0) { backbone[b] = pdb[i]; ++b; } if(strcmp(pdb[i].atom_ptr->atomname,"N")==0) /* new position */ { ++index; k=1; c=1; while(pdb[k].seq_position!=ENDFLAG) { /* find all non-backbone (including CB) atoms for this position */ if(pdb[k].seq_position==pdb[i].seq_position) { if(!(strcmp(pdb[k].atom_ptr->atomname,"N")==0 || strcmp(pdb[k].atom_ptr->atomname,"H")==0 || strcmp(pdb[k].atom_ptr->atomname,"1H")==0 || strcmp(pdb[k].atom_ptr->atomname,"2H")==0 || strcmp(pdb[k].atom_ptr->atomname,"3H")==0 || strcmp(pdb[k].atom_ptr->atomname,"CA")==0 || strcmp(pdb[k].atom_ptr->atomname,"HA")==0 || strcmp(pdb[k].atom_ptr->atomname,"C")==0 || strcmp(pdb[k].atom_ptr->atomname,"O")==0 || strcmp(pdb[k].atom_ptr->atomname,"OT")==0)) { sidechain[index][c] = pdb[k]; ++c; } } ++k; } sidechain[index][c].seq_position=ENDFLAG; sidechain[index][c].atom_ptr=NULL; sidechain[index][c].sasa=ENDFLAG; } ++i; } backbone[b].seq_position=ENDFLAG; backbone[b].atom_ptr=NULL; backbone[b].sasa=ENDFLAG; } /* given num_res residues in mini_pdbATOM pdb, generate the energy_profile_table input_energy_table. input_energy_table == NULL ---> memory will be allocated by this function. */ void generate_energy_profile_table(double ***input_energy_table, mini_pdbATOM *pdb, int num_res) { int i,j, a,b,c; mini_pdbATOM *backbone, **sidechain, *backbone_and_sidechain, *sidechain_and_sidechain; extern int GB_FLAG, SASA_FLAG, COULOMB_FLAG, TORSION_FLAG; SASA_SUM sasa_sum; double **energy_table, dummydouble; ENERGY dummyenergy; /* allocate memory if needed */ backbone = (mini_pdbATOM *)calloc((num_res+2)*6,sizeof(mini_pdbATOM)); backbone_and_sidechain = (mini_pdbATOM *)calloc((num_res+2)*6+MAX_RES_SIZE,sizeof(mini_pdbATOM)); sidechain_and_sidechain = (mini_pdbATOM *)calloc(2*MAX_RES_SIZE+2,sizeof(mini_pdbATOM)); sidechain = (mini_pdbATOM **)calloc(num_res+2,sizeof(mini_pdbATOM *)); for(i=1;i<=num_res;++i) sidechain[i] = (mini_pdbATOM *)calloc(MAX_RES_SIZE,sizeof(mini_pdbATOM)); sidechain[num_res+1]=NULL; if(*input_energy_table==NULL) { energy_table = (double **)calloc(num_res+3,sizeof(double *)); for(i=0;i<=num_res+1;++i) energy_table[i] = (double *)calloc(num_res+3,sizeof(double)); energy_table[num_res+2]=NULL; } else { energy_table = *input_energy_table; } for(i=0;i<=num_res+1;++i) for(j=0;j<=num_res+1;++j) energy_table[i][j]=0.0; /* calculate sasa and born radii for each atom */ energy_calc(pdb,1,1); /* extract sub structures */ extract_sub_structures(pdb, backbone,sidechain); /* calculate internal backbone energy */ dummyenergy = energy_calc(backbone,2,0); energy_table[0][0] = ENERGY_TOTAL_SCALED(dummyenergy); /* calculate sidechain internal energy and sidechain-bkbn energy */ for(i=1;i<=num_res;++i) { /* sidechain internal energy */ dummyenergy = energy_calc(sidechain[i],2,0); energy_table[i][i] = ENERGY_TOTAL_SCALED(dummyenergy); /* make structure of sidechain and backbone */ a=1; b=1; while(backbone[a].seq_position!=sidechain[i][1].seq_position) { backbone_and_sidechain[b] = backbone[a]; ++a; ++b; } while(backbone[a].seq_position==sidechain[i][1].seq_position) { backbone_and_sidechain[b] = backbone[a]; ++a; ++b; } c=1; while(sidechain[i][c].seq_position!=ENDFLAG) { backbone_and_sidechain[b] = sidechain[i][c]; ++c; ++b; } while(backbone[a].seq_position!=ENDFLAG) { backbone_and_sidechain[b] = backbone[a]; ++a; ++b; } backbone_and_sidechain[b].seq_position=ENDFLAG; backbone_and_sidechain[b].atom_ptr=NULL; backbone_and_sidechain[b].sasa=ENDFLAG; dummyenergy = energy_calc(backbone_and_sidechain,2,0); energy_table[i][0] = ENERGY_TOTAL_SCALED(dummyenergy) - (energy_table[0][0] + energy_table[i][i]); energy_table[0][i] = energy_table[i][0]; } /* add backbone sasa-dependent energy to backbone internal energy */ sasa_sum = get_SASA_SUM(backbone); energy_table[0][0] += (double)SASA_FLAG*(GENERAL_ASP*sasa_sum.sasa_total + HYDROPHOBIC_ASP*(sasa_sum.sasa_total - (sasa_sum.N + sasa_sum.O))); /* sidechain pair energies */ for(i=1;i<=num_res;++i) { for(j=(i+1);j<=num_res;++j) { cat_mini_pdbATOM(sidechain[i], sidechain[j], sidechain_and_sidechain); dummyenergy = energy_calc(sidechain_and_sidechain,2,0); energy_table[i][j] = ENERGY_TOTAL_SCALED(dummyenergy) - (energy_table[i][i] + energy_table[j][j]); energy_table[j][i] = energy_table[i][j]; } sasa_sum = get_SASA_SUM(sidechain[i]); energy_table[i][i] += (double)SASA_FLAG*(GENERAL_ASP*sasa_sum.sasa_total + HYDROPHOBIC_ASP*(sasa_sum.sasa_total - (sasa_sum.N + sasa_sum.O))); } /* sum up energies */ energy_table[num_res+1][num_res+1]=0; energy_table[0][num_res+1] = energy_table[0][0]; for(i=1;i<=num_res;++i) { energy_table[i][num_res+1] = energy_table[i][0]; /* side-bkbn */ energy_table[0][num_res+1] += energy_table[0][i]; energy_table[i][num_res+1] += energy_table[i][i]; /* self */ dummydouble=0; for(j=1;j<=num_res;++j) /* side-side pair energies*/ { if(j!=i) dummydouble += energy_table[i][j]; } energy_table[i][num_res+1] += dummydouble; energy_table[num_res+1][i] = energy_table[i][num_res+1]; energy_table[num_res+1][num_res+1] += (energy_table[num_res+1][i] - 0.5*dummydouble); } energy_table[num_res+1][0] = energy_table[0][num_res+1]; energy_table[num_res+1][num_res+1] += energy_table[0][0]; /* backbone internal */ free_memory(backbone); free_memory(backbone_and_sidechain); free_memory(sidechain_and_sidechain); for(i=1;i<=num_res;++i) free_memory(sidechain[i]); free_memory(sidechain); /* set input_energy_table pointer to generated table */ *input_energy_table = energy_table; }