/* EGAD: output_stuff.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 outputting data. */ #include "output_stuff.h" /* prints to file_ptr stream the sum of probabilities of each residuetype in a VARIABLE_POSITION array that has been subjected to mean-field optimization of residue probabilities. composition is an array in which the sum of probabs are listed as ACcDdEeFGHhIKkLMNPQRSTVWYyXx */ void fprintf_composition_vector(FILE *file_ptr, PROTEIN *protein) { int i, i_res, i_res_rot, num_res; double TdS; extern double RT; VARIABLE_POSITION *varPos; varPos = protein->var_pos; i=1; num_res=1; while(varPos[i].seq_position!=ENDFLAG) { while(varPos[i].seq_position!=protein->seqpos_text_map[num_res].seq_position) ++num_res; /* calculate total entropy TdS for this position */ TdS = 0; for(i_res=1;i_res<=varPos[i].number_of_choices;++i_res) { for(i_res_rot=1; i_res_rot<= varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot) if(varPos[i].choice[i_res].lookup_res_ptr->lookupRot[i_res_rot].P[0] >= EPS) TdS += varPos[i].choice[i_res].lookup_res_ptr->lookupRot[i_res_rot].P[0]*log(varPos[i].choice[i_res].lookup_res_ptr->lookupRot[i_res_rot].P[0]); } TdS = -RT*TdS; fprintf(file_ptr, "SCMF_ENTROPY\t%s\tTdS = %f\t", protein->seqpos_text_map[num_res].seqpos_text, TdS); for(i_res=1;i_res<=varPos[i].number_of_choices;++i_res) { fprintf(file_ptr, "%c %f\t", varPos[i].choice[i_res].resparam_ptr->one_letter_code[0], varPos[i].choice[i_res].lookup_res_ptr->P); } fprintf(file_ptr, "\n"); ++i; } } /* output an energy pair table for a given structure to a file. uses energies from the lookup table (JOBTYPE PAIR_ENERGY_TABLE). Use JOBTYPE ENERGY_PROFILE instead, since it uses actual Born radii, sasa, and is faster. This is still kept because this gives the correct energies when local minimization is used */ void output_lookup_table(PROTEIN *protein) { int i,j, a,b,num_res; char *filename, seqpos_string[10]; FILE *output; GENE j_gene, a_gene,b_gene; float E_i_total, temp_energy; LOOKUP_ROTAMER_X_HBOND *lookupRotX_hbond; lookupRotX_hbond = NULL; filename = (char *)calloc(MXLINE_INPUT,sizeof(char)); sprintf(filename, "%s.lookup_table", protein->parameters.output_prefix); output = fopen_file(filename,"w"); /* loads up sidechain structures if they're not already loaded */ protein->final_energy->E_vdw = relaxed_vdw_energy(&protein->final_chr, protein->mini_fixed_atoms, &protein->parameters); fprintf(output,"\tbkbn\t"); protein->final_chr.genes = protein->final_chr.firstgene; num_res=1; while(protein->final_chr.genes->seq_position != ENDFLAG) { while(protein->final_chr.genes->seq_position != protein->seqpos_text_map[num_res].seq_position) ++num_res; sprintf(seqpos_string,"%s",protein->seqpos_text_map[num_res].seqpos_text); fprintf(output,"%s\t",seqpos_string); protein->final_chr.genes = protein->final_chr.genes->nextgene; } protein->final_chr.genes = protein->final_chr.firstgene; fprintf(output,"total\n"); /* bkbn energy */ E_i_total = protein->parameters.fixedatoms_energy; fprintf(output,"bkbn\t%f\t",protein->parameters.fixedatoms_energy); protein->final_chr.genes = protein->final_chr.firstgene; while(protein->final_chr.genes->seq_position != ENDFLAG) { fprintf(output,"%f\t",protein->final_chr.genes->lookupRot->energy_var_fix); E_i_total += protein->final_chr.genes->lookupRot->energy_var_fix; protein->final_chr.genes = protein->final_chr.genes->nextgene; } fprintf(output,"%f\n",E_i_total); i=1; protein->final_chr.genes = protein->final_chr.firstgene; num_res=1; while(protein->final_chr.genes->seq_position != ENDFLAG) { while(protein->final_chr.genes->seq_position != protein->seqpos_text_map[num_res].seq_position) ++num_res; sprintf(seqpos_string,"%s",protein->seqpos_text_map[num_res].seqpos_text); /* sidechain-bkbn plus self energy */ E_i_total = protein->final_chr.genes->lookupRot->energy_var_fix; fprintf(output,"%s\t%f\t",seqpos_string, protein->final_chr.genes->lookupRot->energy_var_fix); j_gene = protein->final_chr.firstgene; j=1; while(j_gene->seq_position != ENDFLAG) { if(j_gene->seq_position != protein->final_chr.genes->seq_position) /* pair energy */ { temp_energy = 0; if(j>i) { a_gene = protein->final_chr.genes; b_gene = j_gene; a = i; b = j; } else { b_gene = protein->final_chr.genes; a_gene = j_gene; b = i; a = j; } temp_energy = rotamer_pair_energy(a_gene, b_gene, &lookupRotX_hbond); if(lookupRotX_hbond!=NULL) free_LOOKUP_ROTAMER_X_HBOND(&lookupRotX_hbond); E_i_total += temp_energy; fprintf(output,"%f\t", temp_energy); } else fprintf(output,"%d\t",ENDFLAG); ++j; j_gene = j_gene->nextgene; } ++i; protein->final_chr.genes = protein->final_chr.genes->nextgene; fprintf(output,"%f\n",E_i_total); } fclose(output); free_memory(filename); } /* prints each energy component */ void fprintf_energy(FILE *output_stream,ENERGY energy) { fprintf(output_stream,"E_vdw = %f\n", energy.E_vdw); fprintf(output_stream,"E_coulomb = %f\n", energy.E_coulomb); fprintf(output_stream,"E_1_4 = %f\n", energy.E_1_4); fprintf(output_stream,"E_born = %f\n", energy.E_born); fprintf(output_stream,"E_pol = %f\n", energy.E_pol); fprintf(output_stream,"E_sasa = %f\n", energy.E_sasa); fprintf(output_stream,"E_hbond = %f\n", energy.E_hbond); fprintf(output_stream,"E_rss = %f\n", energy.E_rss); fprintf(output_stream,"E_specificity = %f\n", energy.E_specificity); fprintf(output_stream,"E_solubility = %f\n", energy.E_solubility); fprintf(output_stream,"E_reference = %f\n", energy.E_reference); fprintf(output_stream,"total = %f\n", ENERGY_TOTAL(energy)); fprintf(output_stream,"E_folded = %f\n", energy.E_folded); fprintf(output_stream,"E_unfolded = %f\n", energy.E_unfolded); fprintf(output_stream,"pseudo_total = %f\n", ENERGY_TOTAL_SCALED(energy)); } void fprintf_residuelevel_sasa(FILE *output_stream, pdbATOM *pdb) { int i, current_seq_position; double side_sasa, bkbn_sasa; i=1; while(pdb[i].seq_position!=ENDFLAG) { current_seq_position = pdb[i].seq_position; side_sasa=0; bkbn_sasa=0; while(current_seq_position == pdb[i].seq_position) { if(pdb[i].atom_ptr->other_info>0 || strcmp(pdb[i].atom_ptr->atomname,"CB")==0 ) side_sasa += pdb[i].sasa; else bkbn_sasa += pdb[i].sasa; ++i; } fprintf(output_stream,"%s\t\t%s\t\t%lf\t%lf\n",pdb[i-1].seqpos_text, pdb[i-1].residuetype, side_sasa,bkbn_sasa+side_sasa); } } void fprintf_sasa(FILE *output_stream,SASA_SUM sasa_sum) { fprintf(output_stream,"sasa_sp3_S: %f\n", sasa_sum.sp3_S); fprintf(output_stream,"sasa_sp2: %f\n", sasa_sum.sp2); fprintf(output_stream,"sasa_O: %f\n", sasa_sum.O); fprintf(output_stream,"sasa_N: %f\n", sasa_sum.N); fprintf(output_stream,"sasa_H: %f\n", sasa_sum.H); fprintf(output_stream,"sasa_hphob: %f\n", sasa_sum.sp2 + sasa_sum.sp3_S + sasa_sum.H); fprintf(output_stream,"sasa_total: %f\n", sasa_sum.sasa_total); fprintf(output_stream,"fraction_sasa_hphob: %f\n", (sasa_sum.sp2 + sasa_sum.sp3_S + sasa_sum.H)/sasa_sum.sasa_total); fprintf(output_stream,"transfer_free_energy_density: %f\n", sasa_sum.E_transfer/sasa_sum.sasa_total); } void fprintf_clashes(FILE *output_stream, pdbATOM *pdb) { int i,j,num_clashes; double r2,r,E; mini_pdbATOM *mini_pdb, *firstAtom, *contactAtom, *twobondAtom, *dihedAtom; VDW_ENERGY vdwE; extern double VDW_CLASH_ENERGY; mini_pdb = (mini_pdbATOM *)calloc(MAX_ATOMS,sizeof(mini_pdbATOM)); pdbATOM_to_mini_pdbATOM(pdb,mini_pdb); fprintf(output_stream,"CLASH\t\tatom1\t\t\t\tatom2\t\t\tr\tE_lennard_jones\tE_vdw_working\n"); i=1; num_clashes=0; while(pdb[i].seq_position!=ENDFLAG) { if(pdb[i].atom_ptr->atom_ptr->atom_class!=ENDFLAG) { j=i+1; while(pdb[j].seq_position!=ENDFLAG) { if(pdb[j].atom_ptr->atom_ptr->atom_class!=ENDFLAG) { r2 = Distance_sqrd(pdb[i].coord, pdb[j].coord); if(get_atom_connectivity(r2, &mini_pdb[i], &mini_pdb[j], &firstAtom, &contactAtom, &twobondAtom, &dihedAtom) == 2) { r = sqrt(r2); VDW_ENERGY_0(vdwE); add_atompair_vdw_energy(&vdwE, mini_pdb[i].atom_ptr, mini_pdb[j].atom_ptr, r2, r); E = sum_VDW_ENERGY(vdwE); if(E > VDW_CLASH_ENERGY) { fprintf(output_stream,"CLASH\t%s\t%s\t%s\t\t%s\t%s\t%s\t\t%lf\t%lf\t%lf\n", pdb[i].seqpos_text, pdb[i].residuetype, pdb[i].atomname, pdb[j].seqpos_text, pdb[j].residuetype, pdb[j].atomname, r, lennard_jones_energy(r2, pdb[i].atom_ptr, pdb[j].atom_ptr), E); ++num_clashes; } } } ++j; } } ++i; } free_memory(mini_pdb); } /* This function outputs all the energies, sequences, and rotamers encoded by the top num_solutions in the CHROMOSOME array to the file name.out */ void output_data(CHROMOSOME *chr, int *num_solutions, const char *name, PROTEIN *protein) { FILE *output; char *filename,*sequence, *var_sequence; int num_sol, n, num_res,p; pdbATOM *pdb; RESPARAM *resparam; resparam = protein->resparam; filename=(char *)calloc(MXLINE_INPUT,sizeof(char)); sequence=(char *)calloc(MAX_RESIDUES,sizeof(char)); var_sequence=(char *)calloc(MAX_RESIDUES,sizeof(char)); sprintf(filename, "%s.out", name); pdb = (pdbATOM *)calloc(MAX_ATOMS, sizeof(pdbATOM)); output = fopen_file(filename,"w"); n=1; num_sol=1; while(num_sol<=(*num_solutions) && chr[n].energy!=ENDFLAG) { fprintf(output,"SOLUTION_# %d\n", num_sol); fprintf(output,"ENERGY\tE_working: %f\n", chr[n].energy); chr[n].genes = chr[n].firstgene; num_res=1; if(chr[n].genes->choice_ptr->lookup_res_ptr!=NULL) { while(chr[n].genes->seq_position!=ENDFLAG) { while(chr[n].genes->seq_position != protein->seqpos_text_map[num_res].seq_position) ++num_res; fprintf(output,"ROTAMER\t %s\t%c\t%s \t ", protein->seqpos_text_map[num_res].seqpos_text,chr[n].genes->varpos_ptr->core_flag, chr[n].genes->choice_ptr->resparam_ptr->one_letter_code); for(p=1;p<=chr[n].genes->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi;++p) { if(chr[n].genes->choice_ptr->resparam_ptr->ligand_flag == 1) fprintf(output,"%.13lf\t ",chr[n].genes->chi[p]); else fprintf(output,"%.1f\t ",chr[n].genes->chi[p]); } fprintf(output,"\n"); chr[n].genes = chr[n].genes->nextgene; } chr[n].genes = chr[n].firstgene; } CHROMOSOME_to_sequence(&chr[n], protein->invar_pos, sequence); fprintf(output,"SEQUENCE_%d\t%s\n",num_sol,sequence); CHROMOSOME_to_variable_sequence(&chr[n], var_sequence); if(var_sequence[0]!=0) fprintf(output,"VARIABLE_SEQUENCE_%d\t%s\n",num_sol,var_sequence); fprintf(output,"END\n"); ++num_sol; ++n; } fclose(output); free_memory(pdb); free_memory(filename); free_memory(sequence); free_memory(var_sequence); } /* This function outputs all the calculated rotamers and the structure encoded by final_chr or final_pdb in a PROTEIN If final_energy == NULL, this function will also calculate the structs and energy, encoded by final_chr */ void output_PROTEIN(PROTEIN *protein, int output_coord_flag) { FILE *output; char dum1[5],dum2[5],dum3[5],dum4[5],dum5[5],remark[6],energy_remark[7],sasa_remark[5],rotamerlabel[8],ligand_dum1[5]; char *filename, *filename2, seq[10], current_chain_id, working_chain_id, *var_sequence; double charge, max_distance, mass, ext_coef; int i, j, n, multichain_flag, nt_flag,ct_flag, num_res, prev_res,q,p, seq_position_eq_1000_flag; div_t chain_id; CARTESIAN centroid; extern int LOCAL_MINIMIZATION_FLAG, GB_FLAG, SASA_FLAG, SOLUBILITY_CUTOFF_FLAG, SPECIFICITY_FLAG,HBOND_SPECIFICITY_FLAG, CHARGE_EQUALITY_FLAG,CHARGES_PH_INDEPENDENT_FLAG, OUTPUT_ENERGY_PER_ATOM_FLAG; extern int INTRA_ROTAMER_FLAG, GET_PID, RESTRAIN_MINIMIZATION_FLAG; extern double TEMPERATURE, IONIC_STRENGTH, PH, INTERNAL_DIELECTRIC,TRANSFER_FREE_ENERGY_DENSITY_CUTOFF; extern double VDW_RADII_SCALE, FRACTION_HYDROPHOBIC_SASA_CUTOFF, OVERALL_CHARGE; extern double HYDROPHOBIC_ASP, GENERAL_ASP, WEIGHT_VDW, WEIGHT_ELECTROSTATICS, WEIGHT_1_4; extern double VDW_LINEAR_START_POINT, VDW_ATTRACTIVE_FACTOR, VDW_REPULSIVE_FACTOR; extern char *INPUTFILENAME; int H,N,C,O,S; ENERGY atom_energy; double E_atom, residue_energy, sidechain_energy, sasa_residue, sasa_side; filename=(char *)calloc(MXLINE_INPUT,sizeof(char)); q=0; multichain_flag=0; if(protein->var_pos !=NULL) { i=1; while(protein->var_pos[i].seq_position!=ENDFLAG) { if(protein->var_pos[i].seq_position > 1000) multichain_flag=1; ++i; } } else { i=1; while(protein->final_pdb[i].seq_position!=ENDFLAG) { if(protein->final_pdb[i].seq_position > 1000) multichain_flag=1; ++i; } } if(protein->chain_id_list!=NULL) multichain_flag=1; strcpy(dum1,"ATOM"); strcpy(ligand_dum1,"ATOM"); strcpy(dum2,"1.00"); strcpy(dum3,"0.00"); strcpy(dum4,""); strcpy(dum5," "); strcpy(remark,"REMARK"); strcpy(energy_remark,"ENERGY"); strcpy(sasa_remark,"SASA"); strcpy(seq,"SEQUENCE"); strcpy(rotamerlabel,"ROTAMER"); sprintf(filename, "%s.pdb.%d", protein->parameters.output_prefix, GET_PID); output = fopen_file(filename,"w"); fprintf(output,"%s\t%s.pdb\n",remark, protein->parameters.output_prefix); fprintf(output,"%s\n", remark); fprintf(output,"TEMPLATE\t%s\n",protein->Template_filename); fprintf(output,"%s\n", remark); fprintf(output,"INPUTFILE\t%s\n",INPUTFILENAME); fprintf(output,"%s\n", remark); fprintf(output,"%s\tTEMPERATURE %f\n", remark, TEMPERATURE); fprintf(output,"%s\tIONIC_STRENGTH %f\n", remark, IONIC_STRENGTH); fprintf(output,"%s\tpH %f\n", remark, PH); fprintf(output,"%s\n", remark); fprintf(output,"%s\tCHARGES_PH_INDEPENDENT_FLAG %d\n", remark, CHARGES_PH_INDEPENDENT_FLAG); fprintf(output,"%s\n", remark); fprintf(output,"%s\tINTERNAL_DIELECTRIC %f\n", remark, INTERNAL_DIELECTRIC); fprintf(output,"%s\tWATER_DIELECTRIC %f\n", remark, WATER_DIELECTRIC); fprintf(output,"%s\n", remark); fprintf(output,"%s\tGENERAL_ASP %f\n", remark, GENERAL_ASP); fprintf(output,"%s\tHYDROPHOBIC_ASP %f\n", remark, HYDROPHOBIC_ASP); fprintf(output,"%s\n", remark); fprintf(output,"%s\tGB_FLAG %d\n", remark, GB_FLAG); fprintf(output,"%s\tSASA_FLAG %d\n", remark, SASA_FLAG); fprintf(output,"%s\n", remark); fprintf(output,"%s\tVDW_RADII_SCALE %f\n", remark, VDW_RADII_SCALE); fprintf(output,"%s\tVDW_ATTRACTIVE_FACTOR %f\n", remark, VDW_ATTRACTIVE_FACTOR); fprintf(output,"%s\tVDW_REPULSIVE_FACTOR %f\n", remark, VDW_REPULSIVE_FACTOR); fprintf(output,"%s\tVDW_LINEAR_START_POINT %f\n", remark, VDW_LINEAR_START_POINT); fprintf(output,"%s\n", remark); fprintf(output,"%s\tWEIGHT_VDW %f\n", remark, WEIGHT_VDW); fprintf(output,"%s\tWEIGHT_ELECTROSTATICS %f\n", remark, WEIGHT_ELECTROSTATICS); fprintf(output,"%s\tWEIGHT_1_4 %f\n", remark, WEIGHT_1_4); fprintf(output,"%s\tOVERALL_ENERGY_SCALE %f\n", remark, OVERALL_ENERGY_SCALE); fprintf(output,"%s\n", remark); fprintf(output,"%s\tINTRA_ROTAMER_FLAG %d\n", remark, INTRA_ROTAMER_FLAG); fprintf(output,"%s\n", remark); if(LOCAL_MINIMIZATION_FLAG!=0 && LOCAL_MINIMIZATION_FLAG!=2) fprintf(output,"%s\tLOCAL_MINIMIZATION_FLAG %d\n", remark, LOCAL_MINIMIZATION_FLAG); if(LOCAL_MINIMIZATION_FLAG==2) fprintf(output,"%s\tMINIMIZE_ROTAMERS_FLAG 1\n", remark); fprintf(output,"%s\n", remark); if(strcmp(protein->parameters.algorithm, "POINT_ENERGY")==0) fprintf(output,"%s\tJOBTYPE POINT_ENERGY\n", remark); else if(strcmp(protein->parameters.algorithm, "MINIMIZE")==0 || strcmp(protein->parameters.algorithm, "POWELL")==0 || strcmp(protein->parameters.algorithm, "DOWNHILL_SIMPLEX")==0 ) { fprintf(output,"%s\tJOBTYPE %s\n", remark, protein->parameters.algorithm); if(protein->parameters.rebuild_backbone_flag == 1) { fprintf(output,"%s\t\tmoving_backbone\t", remark); if(RESTRAIN_MINIMIZATION_FLAG == 0) fprintf(output,"unrestrained\n"); else if(RESTRAIN_MINIMIZATION_FLAG == 1) fprintf(output,"dihedral_restrained %.2lf deg\n",protein->parameters.max_bkbn_delta_dihedral); else if(RESTRAIN_MINIMIZATION_FLAG == 2) fprintf(output,"rmsd_restrained\n"); } } else { if(SOLUBILITY_CUTOFF_FLAG == 1) { fprintf(output,"%s\tFRACTION_HYDROPHOBIC_SASA_CUTOFF %f\n", remark, FRACTION_HYDROPHOBIC_SASA_CUTOFF); fprintf(output,"%s\tTRANSFER_FREE_ENERGY_DENSITY_CUTOFF %f\n", remark, TRANSFER_FREE_ENERGY_DENSITY_CUTOFF); if(CHARGE_EQUALITY_FLAG == 1) { if(OVERALL_CHARGE > 0) fprintf(output,"%s\tOVERALL_CHARGE +%f\n", remark, OVERALL_CHARGE); else fprintf(output,"%s\tOVERALL_CHARGE %f\n", remark, OVERALL_CHARGE); } else if(CHARGE_EQUALITY_FLAG == 0) { fprintf(output,"%s\tOVERALL_CHARGE |%f|\n", remark, OVERALL_CHARGE); } else if(CHARGE_EQUALITY_FLAG == 2) { if(OVERALL_CHARGE > 0) fprintf(output,"%s\tOVERALL_CHARGE >+%f\n", remark, OVERALL_CHARGE); else fprintf(output,"%s\tOVERALL_CHARGE >%f\n", remark, OVERALL_CHARGE); } else if(CHARGE_EQUALITY_FLAG == -2) { if(OVERALL_CHARGE > 0) fprintf(output,"%s\tOVERALL_CHARGE <+%f\n", remark, OVERALL_CHARGE); else fprintf(output,"%s\tOVERALL_CHARGE <%f\n", remark, OVERALL_CHARGE); } else if(CHARGE_EQUALITY_FLAG == 3) fprintf(output,"%s\tOVERALL_CHARGE >|%f|\n", remark, OVERALL_CHARGE); else if(CHARGE_EQUALITY_FLAG == -3) fprintf(output,"%s\tOVERALL_CHARGE <|%f|\n", remark, OVERALL_CHARGE); } fprintf(output,"%s\n", remark); fprintf(output,"%s\tSPECIFICITY_FLAG\t%d\n", remark,SPECIFICITY_FLAG); fprintf(output,"%s\tHBOND_SPECIFICITY_FLAG\t%d\n", remark,HBOND_SPECIFICITY_FLAG); fprintf(output,"%s\n", remark); if(protein->parameters.disk_lookup_table_flag !=0) fprintf(output,"%s\tLOOKUP_TABLE_DIRECTORY %s\n", remark, protein->parameters.lookup_energy_table_directory); fprintf(output,"%s\n", remark); if(strcmp(protein->parameters.algorithm, "GA")==0) { fprintf(output,"%s\tJOBTYPE GENETIC_ALGORITHM\n", remark); fprintf(output,"%s\t\tpopulation %d\n", remark, protein->parameters.pop_size); fprintf(output,"%s\t\tnumber_GA_cycles %d\n", remark, protein->parameters.number_GA_cycles); fprintf(output,"%s\t\tnumber_GA_solutions %d\n", remark, protein->parameters.num_GA_solns_per_cycle); fprintf(output,"%s\t\tconvergence %d\n", remark, protein->parameters.ga_convergence); } else if(strcmp(protein->parameters.algorithm, "SCMF")==0) { fprintf(output,"%s\tJOBTYPE SCMF\n", remark); fprintf_composition_vector(output, protein); fprintf(output,"%s\n", remark); fprintf(output,"%s\tE_scmf: %f\n", energy_remark, protein->E_scmf); } else if(strcmp(protein->parameters.algorithm, "MC")==0) { fprintf(output,"%s\tJOBTYPE MONTE_CARLO\n", remark); fprintf(output,"%s\t\tnumber_MC_cycles %d\n", remark, protein->parameters.number_MC_cycles); fprintf(output,"%s\t\tnumber_MC_solutions %d\n", remark, protein->parameters.number_MC_solns); } else if(strncmp(protein->parameters.algorithm, "MC_GA",5)==0) { fprintf(output,"%s\tJOBTYPE MONTE_CARLO_GENETIC_ALGORITHM\n", remark); fprintf(output,"%s\t\tnumber_MC_cycles %d\n", remark, protein->parameters.number_MC_cycles); fprintf(output,"%s\t\tnumber_MC_solutions %d\n", remark, protein->parameters.number_MC_solns); fprintf(output,"%s\t\tpopulation %d\n", remark, protein->parameters.pop_size); fprintf(output,"%s\t\tnumber_GA_cycles %d\n", remark, protein->parameters.number_GA_cycles); fprintf(output,"%s\t\tnumber_GA_solutions %d\n", remark, protein->parameters.num_GA_solns_per_cycle); fprintf(output,"%s\t\tconvergence %d\n", remark, protein->parameters.ga_convergence); } else if(strcmp(protein->parameters.algorithm, "DEE")==0) { fprintf(output,"%s\tJOBTYPE DEE\n", remark); } else if(strcmp(protein->parameters.algorithm, "FASTER")==0) { fprintf(output,"%s\tJOBTYPE FASTER\n", remark); } else if(strcmp(protein->parameters.algorithm, "ALIGN_STRUCTURES")==0) { fprintf(output,"%s\tJOBTYPE ALIGN_STRUCTURES\n", remark); fprintf(output,"%s\t\trmsd = %lf Angstroms for %d atoms\n",remark,protein->align_rmsd,protein->num_align_points); } } fprintf(output,"%s\n", remark); if(protein->final_chr.genes != NULL) { if(protein->final_pdb == NULL) { protein->final_pdb = (pdbATOM *)calloc(MAX_ATOMS, sizeof(pdbATOM)); CHROMOSOME_to_pdbATOM(&protein->final_chr, protein->fixed_atoms, protein->final_pdb, protein->chain_anchor_bkbn); } } if(protein->final_energy == NULL) final_chr_to_final_pdb_final_energy(protein); pdbATOM_to_sequence(protein->final_pdb, protein->final_sequence, protein->resparam); /* output energies */ fprintf(output,"%s\tE_vdw: %f\n", energy_remark, protein->final_energy->E_vdw); fprintf(output,"%s\tE_coulomb: %f\n", energy_remark, protein->final_energy->E_coulomb); fprintf(output,"%s\tE_1_4: %f\n", energy_remark, protein->final_energy->E_1_4); fprintf(output,"%s\tE_born: %f\n", energy_remark, protein->final_energy->E_born); fprintf(output,"%s\tE_pol: %f\n", energy_remark, protein->final_energy->E_pol); fprintf(output,"%s\tE_sasa: %f\n", energy_remark, protein->final_energy->E_sasa); fprintf(output,"%s\tE_hbond: %f\n", energy_remark, protein->final_energy->E_hbond); fprintf(output,"%s\n",remark); if(fabs(protein->final_energy->TdS) > 1e10) fprintf(output,"%s\tTdS: 0\n",energy_remark); else fprintf(output,"%s\tTdS: %f\n", energy_remark, protein->final_energy->TdS); fprintf(output,"%s\n",remark); fprintf(output,"%s\tE_rss: %f\n", energy_remark, protein->final_energy->E_rss); fprintf(output,"%s\tE_specificity: %f\n", energy_remark, protein->final_energy->E_specificity); fprintf(output,"%s\tE_solubility: %f\n", energy_remark, protein->final_energy->E_solubility); fprintf(output,"%s\n",remark); fprintf(output,"%s\tE_working: %f\n", energy_remark, protein->E_working); fprintf(output,"%s\n",remark); fprintf(output,"%s\tE_structure: %f\n", energy_remark, protein->final_energy->E_structure); fprintf(output,"%s\tE_reference: %f\n", energy_remark, protein->final_energy->E_reference); fprintf(output,"%s\tE_total: %f\n", energy_remark, protein->final_energy->E_total); fprintf(output,"%s\n",remark); fprintf(output,"%s\tE_folded: %f\n", energy_remark, protein->final_energy->E_folded); fprintf(output,"%s\tE_unfolded: %f\n", energy_remark, protein->final_energy->E_unfolded); fprintf(output,"%s\tPseudo_DELTA_G_folding: %f\n", energy_remark, protein->pseudo_dG); if(SASA_FLAG==1) { fprintf(output,"%s\n",remark); fprintf(output,"%s\tsasa_sp3_S: %f\n", sasa_remark, protein->sasa_sum.sp3_S); fprintf(output,"%s\tsasa_sp2: %f\n", sasa_remark, protein->sasa_sum.sp2); fprintf(output,"%s\tsasa_O: %f\n", sasa_remark, protein->sasa_sum.O); fprintf(output,"%s\tsasa_N: %f\n", sasa_remark, protein->sasa_sum.N); fprintf(output,"%s\tsasa_H: %f\n", sasa_remark, protein->sasa_sum.H); fprintf(output,"%s\tsasa_hphob: %f\n", sasa_remark, protein->sasa_sum.sp2 + protein->sasa_sum.sp3_S + protein->sasa_sum.H); fprintf(output,"%s\tsasa_total: %f\n", sasa_remark, protein->sasa_sum.sasa_total); fprintf(output,"%s\tfraction_sasa_hphob: %f\n", sasa_remark, (protein->sasa_sum.sp2 + protein->sasa_sum.sp3_S + protein->sasa_sum.H)/protein->sasa_sum.sasa_total); fprintf(output,"%s\ttransfer_free_energy_density: %f\n", sasa_remark, protein->sasa_sum.E_transfer/protein->sasa_sum.sasa_total); } fprintf(output,"%s\n",remark); if(protein->transform_matrix!=NULL) { fprintf(output,"TRANSFORM_MATRIX\t%.13lf %.13lf %.13lf %.13lf\n", protein->transform_matrix[1],protein->transform_matrix[2],protein->transform_matrix[3],protein->transform_matrix[4]); fprintf(output,"TRANSFORM_MATRIX\t%.13lf %.13lf %.13lf %.13lf\n", protein->transform_matrix[5],protein->transform_matrix[6],protein->transform_matrix[7],protein->transform_matrix[8]); fprintf(output,"TRANSFORM_MATRIX\t%.13lf %.13lf %.13lf %.13lf\n", protein->transform_matrix[9],protein->transform_matrix[10],protein->transform_matrix[11],protein->transform_matrix[12]); fprintf(output,"%s\n",remark); } charge = 0; /* output rotamers */ if(protein->final_chr.genes != NULL) { fprintf(output,"%s\tSEARCH_SPACE: 10^%lf sequences\t10^%lf rotamer_combinations\n",remark, protein->parameters.log10_seq_combinations, protein->parameters.log10_rotamer_combinations); fprintf(output,"%s\n",remark); protein->final_chr.genes = protein->final_chr.firstgene; i=1; num_res = 1; while(protein->final_chr.genes->seq_position!=ENDFLAG) { while(protein->final_chr.genes->seq_position != protein->seqpos_text_map[num_res].seq_position) ++num_res; if(protein->final_chr.genes->varpos_ptr->core_flag == 'c' || protein->final_chr.genes->varpos_ptr->core_flag == 's' || protein->final_chr.genes->varpos_ptr->core_flag == 'i' || protein->final_chr.genes->varpos_ptr->core_flag == 'l') { if(protein->final_chr.genes->varpos_ptr->fixed_flag == 1) fprintf(output,"%s\t %s\t%c.fix\t%s \t ",rotamerlabel, protein->seqpos_text_map[num_res].seqpos_text, protein->final_chr.genes->varpos_ptr->core_flag, protein->final_chr.genes->choice_ptr->resparam_ptr->one_letter_code); else { if(protein->final_chr.genes->varpos_ptr->number_of_choices > 1) fprintf(output,"%s\t %s\t%c.des\t%s \t ",rotamerlabel, protein->seqpos_text_map[num_res].seqpos_text, protein->final_chr.genes->varpos_ptr->core_flag, protein->final_chr.genes->choice_ptr->resparam_ptr->one_letter_code); else fprintf(output,"%s\t %s\t%c.rot\t%s \t ",rotamerlabel, protein->seqpos_text_map[num_res].seqpos_text, protein->final_chr.genes->varpos_ptr->core_flag, protein->final_chr.genes->choice_ptr->resparam_ptr->one_letter_code); } for(p=1;p<=protein->final_chr.genes->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi;++p) { if(protein->final_chr.genes->choice_ptr->resparam_ptr->ligand_flag==1) fprintf(output,"%.13lf\t ",protein->final_chr.genes->chi[p]); else fprintf(output,"%.1f\t ",protein->final_chr.genes->chi[p]); } fprintf(output,"\n"); } ++i; protein->final_chr.genes = protein->final_chr.genes->nextgene; } protein->final_chr.genes = protein->final_chr.firstgene; } /* fprintf(output,"%s\n",remark); fprintf(output,"CHARGE\t%lf\n",charge); */ i=1; n=1; while(protein->final_pdb[i].seq_position!=ENDFLAG) { if(strcmp(protein->final_pdb[i].residuetype, "NTE")==0) { n=i; while(protein->final_pdb[i].seq_position%1000==1) { if(protein->final_pdb[i].atom_ptr->atomnumber == 8) { q=i; } ++i; } j=i-1; for(i=n;i<=j;++i) strcpy(protein->final_pdb[i].residuetype, protein->final_pdb[q].residuetype); i=j+1; } ++i; } fprintf(output,"%s\n",remark); i=1; max_distance=-10; while(strcmp(protein->final_pdb[i].residuetype,"END")!=0) { j=i+1; while(strcmp(protein->final_pdb[j].residuetype,"END")!=0) { if(max_distance <= Distance_sqrd(protein->final_pdb[i].coord,protein->final_pdb[j].coord)) { max_distance = Distance_sqrd(protein->final_pdb[i].coord,protein->final_pdb[j].coord); } ++j; } ++i; } /* print out the sequence(s) */ if(protein->final_chr.firstgene!=NULL) { var_sequence = (char *)calloc(MAXLINE,sizeof(char)); CHROMOSOME_to_variable_sequence(&protein->final_chr, var_sequence); if(var_sequence[0]!=0) fprintf(output,"VARIABLE_SEQUENCE\t%s\n",var_sequence); free_memory(var_sequence); } i=1; seq_position_eq_1000_flag=0; if(multichain_flag != 0) { while(strcmp(protein->final_pdb[i].residuetype,"END")!=0) { chain_id = div(protein->final_pdb[i].seq_position, 1000); if(chain_id.rem == 0) seq_position_eq_1000_flag=1; if(seq_position_eq_1000_flag==1) { chain_id.quot = chain_id.quot - 1; chain_id.rem = chain_id.rem + 1000; } current_chain_id = protein->chain_id_list[chain_id.quot+1]; working_chain_id = current_chain_id; fprintf(output,"%s %c\t", seq, current_chain_id); while(strcmp(protein->final_pdb[i].residuetype,"END")!=0 && working_chain_id == current_chain_id) { n=1; while(strcmp(protein->final_pdb[i].residuetype, protein->resparam[n].residuetype)!=0) ++n; if(strcmp(protein->final_pdb[i].residuetype, "CTE")!=0 && strcmp(protein->final_pdb[i].residuetype, "NTE")!=0 ) fprintf(output,"%c", protein->resparam[n].one_letter_code[0]); else { while(strcmp(protein->final_pdb[i].residuetype, "CTE")==0 || strcmp(protein->final_pdb[i].residuetype, "NTE")==0) ++i; n=1; while(strcmp(protein->final_pdb[i].residuetype, protein->resparam[n].residuetype)!=0) ++n; fprintf(output,"%c", protein->resparam[n].one_letter_code[0]); } j=i; /* advance to next position */ while(protein->final_pdb[i].seq_position==protein->final_pdb[j].seq_position) ++i; if(strcmp(protein->final_pdb[i].residuetype,"END")!=0) { chain_id = div(protein->final_pdb[i].seq_position, 1000); if(chain_id.rem == 0) seq_position_eq_1000_flag=1; if(seq_position_eq_1000_flag==1) { chain_id.quot = chain_id.quot - 1; chain_id.rem = chain_id.rem + 1000; } current_chain_id = protein->chain_id_list[chain_id.quot+1]; if(working_chain_id == current_chain_id) /* same chain */ if((protein->final_pdb[i].seq_position - protein->final_pdb[i-1].seq_position) > 1) /* missing residues */ { n=chain_id.rem; /* current seq pos in this chain */ chain_id = div(protein->final_pdb[i-1].seq_position, 1000); if(chain_id.rem == 0) seq_position_eq_1000_flag=1; if(seq_position_eq_1000_flag==1) { chain_id.quot = chain_id.quot - 1; chain_id.rem = chain_id.rem + 1000; } j=chain_id.rem; /* seq pos of prev. residue in this chain */ n=n-j; for(j=1;jfinal_pdb[i].residuetype,"END")!=0) { n=1; while(strcmp(protein->final_pdb[i].residuetype, protein->resparam[n].residuetype)!=0) ++n; if(strcmp(protein->final_pdb[i].residuetype, "CTE")!=0 && strcmp(protein->final_pdb[i].residuetype, "NTE")!=0 ) { fprintf(output,"%c", protein->resparam[n].one_letter_code[0]); charge += protein->resparam[n].overall_charge_pH_avg; if(protein->resparam[n].one_letter_code[0] == 'W') ext_coef += 5690.0; else if(protein->resparam[n].one_letter_code[0] == 'Y' || protein->resparam[n].one_letter_code[0] == 'y') ext_coef += 1280.0; } else { if(strcmp(protein->final_pdb[i].residuetype, "NTE")==0) nt_flag=1; if(strcmp(protein->final_pdb[i].residuetype, "CTE")==0) ct_flag=1; /* bkbn masses */ while(strcmp(protein->final_pdb[i].residuetype, "CTE")==0 || strcmp(protein->final_pdb[i].residuetype, "NTE")==0) { if(protein->final_pdb[i].atom_ptr->atom_ptr->torsiontype[0] == 'H' && protein->final_pdb[i].atom_ptr->atom_ptr->torsiontype[1] != 'U') mass += 1.0079; else if(protein->final_pdb[i].atom_ptr->atom_ptr->torsiontype[0] == 'C') mass += 12.011; else if(protein->final_pdb[i].atom_ptr->atom_ptr->torsiontype[0] == 'N') mass += 14.00674; else if(protein->final_pdb[i].atom_ptr->atom_ptr->torsiontype[0] == 'O') mass += 15.9994; else if(protein->final_pdb[i].atom_ptr->atom_ptr->torsiontype[0] == 'S') mass += 32.06; ++i; } n=1; while(strcmp(protein->final_pdb[i].residuetype, protein->resparam[n].residuetype)!=0) ++n; fprintf(output,"%c", protein->resparam[n].one_letter_code[0]); charge += protein->resparam[n].overall_charge_pH_avg; if(protein->resparam[n].one_letter_code[0] == 'W') ext_coef += 5690.0; else if(protein->resparam[n].one_letter_code[0] == 'Y' || protein->resparam[n].one_letter_code[0] == 'y') ext_coef += 1280.0; } j=i; n=0; while(protein->final_pdb[i].seq_position==protein->final_pdb[j].seq_position) { if(protein->final_pdb[i].atom_ptr->atom_ptr->torsiontype[0] == 'H' && protein->final_pdb[i].atom_ptr->atom_ptr->torsiontype[1] != 'U') { mass += 1.0079; ++H; } else if(protein->final_pdb[i].atom_ptr->atom_ptr->torsiontype[0] == 'C') { mass += 12.011; ++C; } else if(protein->final_pdb[i].atom_ptr->atom_ptr->torsiontype[0] == 'N') { mass += 14.00674; ++N; } else if(protein->final_pdb[i].atom_ptr->atom_ptr->torsiontype[0] == 'O') { mass += 15.9994; ++O; if(strcmp(protein->final_pdb[i].atomname,"O")==0) ++n; } else if(protein->final_pdb[i].atom_ptr->atom_ptr->torsiontype[0] == 'S') { mass += 32.06; ++S; } ++i; } if(n>1) /* two O for a residue; presumably at the C-term; for when a C-term oxygen is present but CTE_FLAG=0 */ ct_flag=1; } fprintf(output,"\n%s\n",remark); if(nt_flag == 0) mass += 2*1.0079; if(ct_flag == 0) mass += 15.9994; mass = mass-charge; /* to equal the expasy values for ionizable residues */ fprintf(output,"MASS %.2f\t\tw/ M: %.2f\tw/ GA: %.2f\tw/ MGA: %.2f\n",mass,mass+131.2,mass+128.1,mass+259.3); fprintf(output,"EXTINCTION_COEFFICIENT eps_280nm: %.1f\n",ext_coef); fprintf(output,"CHARGE at pH %.2lf: %.2f\n",PH,charge); } if(protein->super_chain_id_list!=NULL) /* print SUPER_CHAIN identifiers */ { i=1; while(protein->super_chain_id_list[i].chain_id!=NULL) { fprintf(output,"SUPER_CHAIN\t"); n=1; while(protein->super_chain_id_list[i].chain_id[n]!='\0') { fprintf(output,"%c ",protein->super_chain_id_list[i].chain_id[n] ); ++n; } fprintf(output,"\n"); ++i; } fprintf(output,"%s\n",remark); } if(strcmp(protein->parameters.algorithm, "MC_PAIR_ENERGY_TABLE")==0) { output_lookup_table(protein); } if(strcmp(protein->parameters.algorithm, "MC_ENERGY_PROFILE_TABLE")==0) { generate_energy_profile_table(&protein->energy_profile_table, protein->final_minipdb, protein->num_res); output_energy_profile_table(protein->energy_profile_table, protein->parameters.output_prefix, protein->num_res, protein->seqpos_text_map,protein->final_sequence); } fprintf(output,"%s\n",remark); fprintf(output,"%s\tnumber_of_hbonds: %d\t%d\t%d\t%d\t\ttotal side-side side-bkbn bkbn-bkbn\n", "HBOND", protein->final_energy->num_hbond, protein->final_energy->num_ss_hbond, protein->final_energy->num_sb_hbond, protein->final_energy->num_bb_hbond ); fprintf(output,"%s\tnumber_of_unsatisfied_hbond_groups:\t%d\t%d\t%d\t\ttotal side bkbn\n","HBOND",protein->final_energy->num_unsatisfied_hbond,protein->final_energy->num_side_unsatisfied_hbond,protein->final_energy->num_bkbn_unsatisfied_hbond); fprintf(output,"%s\n",remark); fprintf(output,"%s\tnumber_of_clashes:\t%d\n","CLASH",protein->final_energy->clashes); /* if output_coord_flag == 0, close the file now; saves time for processes where we don't care about the structure */ if(output_coord_flag == 1) { fprintf_clashes(output,protein->final_pdb); centroid = get_centroid(protein->final_minipdb); fprintf(output,"%s\n",remark); fprintf(output,"LENGTH\t%f\n",sqrt(max_distance)); fprintf(output,"CENTROID\t%.3f\t%7.3f\t%7.3f\n",centroid.x, centroid.y, centroid.z); fprintf(output,"%s\n",remark); fprintf(output,"%s\n",remark); fprintf(output,"%s occupancy = SASA B-factor = born radius E_atom '!' = unsatisfied hbond group @ E_sidechain E_residue SASA_sidechain SASA_residue\n",remark); fprintf(output,"%s\n",remark); /* print out the pdbATOM array */ /* rename C-terminal oxygens OXT */ i=1; E_atom=0; residue_energy=0; sidechain_energy=0; sasa_residue=0; sasa_side=0; while(protein->final_pdb[i].seq_position!=ENDFLAG) { n=protein->final_pdb[i].seq_position; j=0; while(protein->final_pdb[i].seq_position==n) { if(strcmp(protein->final_pdb[i].atomname,"O")==0) ++j; ++i; } if(j==2) /* two atoms named O for this residue; rename last one to OXT */ { --i; while(strcmp(protein->final_pdb[i].atomname,"O")!=0) --i; strcpy(protein->final_pdb[i].atomname,"OXT"); while(protein->final_pdb[i].seq_position==n) ++i; } } i=1; n=0; num_res=0; prev_res=0; seq_position_eq_1000_flag=0; while(strcmp(protein->final_pdb[i].residuetype,"END")!=0) { if(protein->final_pdb[i].seq_position != prev_res) { prev_res = protein->final_pdb[i].seq_position; ++num_res; } if(strstr(protein->final_pdb[i].atom_ptr->atom_ptr->atomtype, "U")==0) { if(strcmp(protein->final_pdb[i].residuetype, "HIN")==0) strcpy(protein->final_pdb[i].residuetype, "HIS"); if(strcmp(protein->final_pdb[i].residuetype, "CYD")==0) strcpy(protein->final_pdb[i].residuetype, "CYS"); ++n; if(strcmp(protein->final_pdb[i].residuetype, "CTE")==0) { j=i; while(strcmp(protein->final_pdb[j].residuetype, "CTE")==0) { ++j; if(strcmp(protein->final_pdb[j].atomname,"2HA")==0) strcpy(protein->final_pdb[j].residuetype,"GLY"); } strcpy(protein->final_pdb[i].residuetype,protein->final_pdb[j].residuetype); } p=1; while(strcmp(protein->final_pdb[i].residuetype, protein->resparam[p].residuetype)!=0) ++p; if(protein->resparam[p].ligand_flag==0) fprintf(output,"%4s %6d",dum1,n); else { fprintf(output,"%4s %6d",ligand_dum1,n); chain_id = div(protein->final_pdb[i].seq_position, 1000); protein->chain_id_list[chain_id.quot+1] = ' '; } if(multichain_flag==0) { if(protein->final_pdb[i].seq_position<1000) { if(strncmp(protein->final_pdb[i].atomname,"9",1)<=0 || strncmp(protein->final_pdb[i].atomname,"HH",2)==0) fprintf(output,"%s %-4s %-3s %5d %11.3f %7.3f %7.3f %5.2f %.2f", dum4,protein->final_pdb[i].atomname,protein->final_pdb[i].residuetype,protein->final_pdb[i].seq_position, protein->final_pdb[i].coord.x,protein->final_pdb[i].coord.y,protein->final_pdb[i].coord.z, protein->final_pdb[i].sasa, protein->final_pdb[i].born_radius); else fprintf(output,"%s %-3s %-3s %5d %11.3f %7.3f %7.3f %5.2f %.2f", dum5,protein->final_pdb[i].atomname,protein->final_pdb[i].residuetype,protein->final_pdb[i].seq_position, protein->final_pdb[i].coord.x,protein->final_pdb[i].coord.y,protein->final_pdb[i].coord.z, protein->final_pdb[i].sasa, protein->final_pdb[i].born_radius); } else { if(strncmp(protein->final_pdb[i].atomname,"9",1)<=0 || strncmp(protein->final_pdb[i].atomname,"HH",2)==0) fprintf(output,"%s %-4s %-3s%5d %11.3f %7.3f %7.3f %5.2f %.2f", dum4,protein->final_pdb[i].atomname,protein->final_pdb[i].residuetype,protein->final_pdb[i].seq_position, protein->final_pdb[i].coord.x,protein->final_pdb[i].coord.y,protein->final_pdb[i].coord.z, protein->final_pdb[i].sasa, protein->final_pdb[i].born_radius); else fprintf(output,"%s %-3s %-3s%5d %11.3f %7.3f %7.3f %5.2f %.2f", dum5,protein->final_pdb[i].atomname,protein->final_pdb[i].residuetype,protein->final_pdb[i].seq_position, protein->final_pdb[i].coord.x,protein->final_pdb[i].coord.y,protein->final_pdb[i].coord.z, protein->final_pdb[i].sasa, protein->final_pdb[i].born_radius); } } else { chain_id = div(protein->final_pdb[i].seq_position, 1000); if(chain_id.rem == 0) seq_position_eq_1000_flag=1; if(seq_position_eq_1000_flag==1) { chain_id.quot = chain_id.quot - 1; chain_id.rem = chain_id.rem + 1000; } if(protein->wacky_numbering_list == NULL) { if(chain_id.rem < 1000) { if(strncmp(protein->final_pdb[i].atomname,"9",1)<=0 || strncmp(protein->final_pdb[i].atomname,"HH",2)==0) fprintf(output,"%s %-4s %-3s %c %3d %11.3f %7.3f %7.3f %5.2f %.2f", dum4,protein->final_pdb[i].atomname,protein->final_pdb[i].residuetype,protein->chain_id_list[chain_id.quot+1], chain_id.rem, protein->final_pdb[i].coord.x,protein->final_pdb[i].coord.y,protein->final_pdb[i].coord.z, protein->final_pdb[i].sasa, protein->final_pdb[i].born_radius); else fprintf(output,"%s %-3s %-3s %c %3d %11.3f %7.3f %7.3f %5.2f %.2f", dum5,protein->final_pdb[i].atomname,protein->final_pdb[i].residuetype,protein->chain_id_list[chain_id.quot+1], chain_id.rem, protein->final_pdb[i].coord.x,protein->final_pdb[i].coord.y,protein->final_pdb[i].coord.z, protein->final_pdb[i].sasa, protein->final_pdb[i].born_radius); } else { if(strncmp(protein->final_pdb[i].atomname,"9",1)<=0 || strncmp(protein->final_pdb[i].atomname,"HH",2)==0) fprintf(output,"%s %-4s %-3s %c%3d %11.3f %7.3f %7.3f %5.2f %.2f", dum4,protein->final_pdb[i].atomname,protein->final_pdb[i].residuetype,protein->chain_id_list[chain_id.quot+1], chain_id.rem, protein->final_pdb[i].coord.x,protein->final_pdb[i].coord.y,protein->final_pdb[i].coord.z, protein->final_pdb[i].sasa, protein->final_pdb[i].born_radius); else fprintf(output,"%s %-3s %-3s %c%3d %11.3f %7.3f %7.3f %5.2f %.2f", dum5,protein->final_pdb[i].atomname,protein->final_pdb[i].residuetype,protein->chain_id_list[chain_id.quot+1], chain_id.rem, protein->final_pdb[i].coord.x,protein->final_pdb[i].coord.y,protein->final_pdb[i].coord.z, protein->final_pdb[i].sasa, protein->final_pdb[i].born_radius); } } else { if(protein->wacky_numbering_list[num_res] != ' ') { if(strncmp(protein->final_pdb[i].atomname,"9",1)<=0 || strncmp(protein->final_pdb[i].atomname,"HH",2)==0) fprintf(output,"%s %-4s %-3s %c %3d%c %10.3f %7.3f %7.3f %5.2f %.2f", dum4,protein->final_pdb[i].atomname,protein->final_pdb[i].residuetype, protein->chain_id_list[chain_id.quot+1], chain_id.rem - protein->wacky_numbering_list[num_res], protein->wacky_numbering_list[num_res], protein->final_pdb[i].coord.x,protein->final_pdb[i].coord.y,protein->final_pdb[i].coord.z, protein->final_pdb[i].sasa, protein->final_pdb[i].born_radius); else fprintf(output,"%s %-3s %-3s %c %3d%c %10.3f %7.3f %7.3f %5.2f %.2f", dum5,protein->final_pdb[i].atomname,protein->final_pdb[i].residuetype,protein->chain_id_list[chain_id.quot+1], chain_id.rem - protein->wacky_numbering_list[num_res], protein->wacky_numbering_list[num_res], protein->final_pdb[i].coord.x,protein->final_pdb[i].coord.y,protein->final_pdb[i].coord.z, protein->final_pdb[i].sasa, protein->final_pdb[i].born_radius); } else { if(strncmp(protein->final_pdb[i].atomname,"9",1)<=0 || strncmp(protein->final_pdb[i].atomname,"HH",2)==0) fprintf(output,"%s %-4s %-3s %c %3d %11.3f %7.3f %7.3f %5.2f %.2f", dum4,protein->final_pdb[i].atomname,protein->final_pdb[i].residuetype,protein->chain_id_list[chain_id.quot+1], chain_id.rem, protein->final_pdb[i].coord.x,protein->final_pdb[i].coord.y,protein->final_pdb[i].coord.z, protein->final_pdb[i].sasa, protein->final_pdb[i].born_radius); else fprintf(output,"%s %-3s %-3s %c %3d %11.3f %7.3f %7.3f %5.2f %.2f", dum5,protein->final_pdb[i].atomname,protein->final_pdb[i].residuetype,protein->chain_id_list[chain_id.quot+1], chain_id.rem, protein->final_pdb[i].coord.x,protein->final_pdb[i].coord.y,protein->final_pdb[i].coord.z, protein->final_pdb[i].sasa,protein->final_pdb[i].born_radius ); } } } if(protein->resparam[p].ligand_flag==1) /* reset */ protein->chain_id_list[chain_id.quot+1] = 'l'; if(OUTPUT_ENERGY_PER_ATOM_FLAG == 1) atom_energy = energy_for_atom_of_interest(protein->final_minipdb, i, GB_FLAG, SASA_FLAG, TORSION_FLAG); else ENERGY_0(atom_energy); E_atom = ENERGY_TOTAL(atom_energy); fprintf(output,"\t%lf", E_atom); if(protein->final_pdb[i].hbond_satisfied == 0) fprintf(output," !"); residue_energy += E_atom; sasa_residue += protein->final_pdb[i].sasa; if(protein->final_pdb[i].atom_ptr->other_info>0 || strcmp(protein->final_pdb[i].atom_ptr->atomname,"CB")==0) { sidechain_energy += E_atom; sasa_side += protein->final_pdb[i].sasa; } if(protein->final_pdb[i].seq_position!=protein->final_pdb[i+1].seq_position) { if(protein->final_pdb[i].hbond_satisfied == 0) fprintf(output,"@ %lf %lf %lf %lf",sidechain_energy, residue_energy, sasa_side, sasa_residue); else fprintf(output," @ %lf %lf %lf %lf",sidechain_energy, residue_energy, sasa_side, sasa_residue); residue_energy = 0; sidechain_energy = 0; sasa_side=0; sasa_residue=0; } else /* the last atom for GLU and ASP is an HU virtual atom */ { if(protein->final_pdb[i].seq_position!=protein->final_pdb[i+2].seq_position) { if(strstr(protein->final_pdb[i+1].atom_ptr->atomtype,"U")!=0) { if(OUTPUT_ENERGY_PER_ATOM_FLAG == 1) atom_energy = energy_for_atom_of_interest(protein->final_minipdb, i+1, GB_FLAG, SASA_FLAG, TORSION_FLAG); else ENERGY_0(atom_energy); E_atom = ENERGY_TOTAL(atom_energy); residue_energy += E_atom; sidechain_energy += E_atom; if(protein->final_pdb[i].hbond_satisfied == 0) fprintf(output,"@ %lf %lf %lf %lf",sidechain_energy, residue_energy, sasa_side, sasa_residue); else fprintf(output," @ %lf %lf %lf %lf",sidechain_energy, residue_energy, sasa_side, sasa_residue); residue_energy = 0; sidechain_energy = 0; sasa_side=0; sasa_residue=0; } } } fprintf(output,"\n"); } ++i; } } fprintf(output,"END\n"); filename2 = (char *)calloc(MAXLINE,sizeof(char)); sprintf(filename2, "%s.pdb", protein->parameters.output_prefix); mv_file(filename, filename2); free_memory(filename); free_memory(filename2); fclose(output); } /* prints coords in pdb to standard out; for debugging */ void printf_minipdb(mini_pdbATOM pdb[]) { int i; char dum1[5],dum2[5],dum3[5]; strcpy(dum1,"ATOM"); strcpy(dum2,"1.00"); strcpy(dum3,"0.00"); i=1; while(pdb[i].seq_position!=ENDFLAG) { if(strncmp(pdb[i].atom_ptr->atomname,"9",1)<=0 || strncmp(pdb[i].atom_ptr->atomname,"HH",2)==0) printf("%-4s %s %5d %11.3f %7.3f %7.3f %f %f\n", pdb[i].atom_ptr->atomname,pdb[i].atom_ptr->atomtype,pdb[i].seq_position, pdb[i].coord.x,pdb[i].coord.y,pdb[i].coord.z, pdb[i].born_radius, pdb[i].sasa ); else printf("%-3s %s %5d %11.3f %7.3f %7.3f %f %f\n", pdb[i].atom_ptr->atomname,pdb[i].atom_ptr->atomtype,pdb[i].seq_position, pdb[i].coord.x,pdb[i].coord.y,pdb[i].coord.z, pdb[i].born_radius, pdb[i].sasa ); ++i; } } /* This function prints a pdbATOM array to the screen in the pdb format; for debugging */ void fprintf_pdb(FILE *output_stream, pdbATOM pdb[]) { int i; char dum1[5],dum2[5],dum3[5]; strcpy(dum1,"ATOM"); strcpy(dum2,"1.00"); strcpy(dum3,"0.00"); i=1; while(strcmp(pdb[i].residuetype,"END")!=0) { fprintf(output_stream,"%4s %6d",dum1,i); if(strncmp(pdb[i].atomname,"9",1)<=0 || strncmp(pdb[i].atomname,"HH",2)==0) fprintf(output_stream,"%s %-4s %-3s %5d %11.3f %7.3f %7.3f %6.3f %.3f\n", "",pdb[i].atomname, pdb[i].residuetype,pdb[i].seq_position, pdb[i].coord.x,pdb[i].coord.y,pdb[i].coord.z, 1.0, 0.0 /* pdb[i].born_radius, pdb[i].sasa */ ); else fprintf(output_stream,"%s %-3s %-3s %5d %11.3f %7.3f %7.3f %6.3f %.3f\n", " ", pdb[i].atomname, pdb[i].residuetype,pdb[i].seq_position, pdb[i].coord.x,pdb[i].coord.y,pdb[i].coord.z, 1.0, 0.0 /* pdb[i].born_radius, pdb[i].sasa */ ); ++i; } } /* This function inputs a CHROMOSOME and the fixed backbone atoms It returns the vdw energy based on relaxing sidechains via rotamerlets and methyl-lets This requires that all the sidechain coords be pre-calculated by generate_lookup_table and linked appropriately */ double relaxed_vdw_energy(CHROMOSOME *chr, mini_pdbATOM *fixed_atoms, PARAMETERS *parameters) { int i,j, h, k, x, q, y; GENE i_gene; /* marker for position i on the CHROMOSOME */ double vdw_energy, dummy_vdw_energy1, energy_vdw; int best_rotamerlet_i, best_rotamerlet_j; double MeH_MeH_energy[MAX_NUM_METHYL][4][MAX_NUM_METHYL][4], nonMeH_nonMeH_energy; double MeH_nonMeH_energy_i[MAX_NUM_METHYL][4], MeH_nonMeH_energy_j[MAX_NUM_METHYL][4]; ROTAMERLET *rotamerlet_Me, *rotamerlet_non_Me; RESPARAM *resparam_ptr_Me, *resparam_ptr_non_Me; int conf_i, conf_j; double total_vdw_energy; var_var_ENERGY var_var_ENERGy; var_fix_vdw_ENERGY dummy_var_fix_vdw_energy; ENERGY energy; var_fix_ENERGY var_fix_ENERGy; int best_methyl, best_rotamerlet; mini_pdbATOM *bkbn_atoms, *sideAtoms; extern int **BASE_THREE_CONVERT; char *structure_filename; bool do_these_hbond_flag, *i_hbond_status, *j_hbond_status; extern int LOCAL_MINIMIZATION_FLAG; /* MINIMIZE_METHYL_FLAG; */ best_rotamerlet_i = ENDFLAG; best_rotamerlet_j = ENDFLAG; best_rotamerlet = ENDFLAG; sideAtoms = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM)); bkbn_atoms = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM)); structure_filename = (char *)calloc(MXLINE_INPUT, sizeof(char)); i_hbond_status = (bool *)calloc(MAX_RES_SIZE,sizeof(bool)); j_hbond_status = (bool *)calloc(MAX_RES_SIZE,sizeof(bool)); sprintf(structure_filename,"coordinates"); create_directory(structure_filename, parameters->lookup_energy_table_directory); i=1; j=1; chr->genes = chr->firstgene; while(chr->genes->seq_position!=ENDFLAG) { while(fixed_atoms[i].seq_position!=chr->genes->seq_position && fixed_atoms[i].atom_ptr!=NULL) { bkbn_atoms[j] = fixed_atoms[i]; ++i; ++j; } while(fixed_atoms[i].seq_position == chr->genes->seq_position && fixed_atoms[i].atom_ptr!=NULL) { if(fixed_atoms[i].atom_ptr->atomnumber!=8 && fixed_atoms[i].atom_ptr->other_info<0) { bkbn_atoms[j] = fixed_atoms[i]; ++j; } ++i; } chr->genes = chr->genes->nextgene; } while(fixed_atoms[i].atom_ptr!=NULL) { bkbn_atoms[j] = fixed_atoms[i]; ++j; ++i; } bkbn_atoms[j] = fixed_atoms[i]; energy = energy_calc(bkbn_atoms, 0, 0); total_vdw_energy = energy.E_vdw; chr->genes = chr->firstgene; i=1; j=1; while(chr->genes->seq_position!=ENDFLAG) { i_gene = chr->genes; total_vdw_energy += iROT->rotamer->energy.E_vdw; if(LOCAL_MINIMIZATION_FLAG == 1) { var_var_ENERGy.E_vdw = DBL_MAX; for(h=1;h<=iRES->rotamerlib_ptr->num_rotamerlets;++h) { if(i_gene->choice_ptr->lookup_res_ptr->lookupRot[1].sideAtoms == NULL) { if(i_gene->choice_ptr->lookup_res_ptr->fixed_flag == 1) { sprintf(structure_filename, "%s/coordinates/%d/%s.%d.structure.fixed", parameters->lookup_energy_table_directory, i_gene->seq_position, i_gene->choice_ptr->resparam_ptr->residuetype, i_gene->seq_position); } else { sprintf(structure_filename, "%s/coordinates/%d/%s.%d.structure", parameters->lookup_energy_table_directory, i_gene->seq_position, i_gene->choice_ptr->resparam_ptr->residuetype, i_gene->seq_position); } load_coordinates_lookupRes(i_gene->choice_ptr->lookup_res_ptr, i_gene->choice_ptr->resparam_ptr, structure_filename); } dummy_var_fix_vdw_energy = var_fixed_vdw_screen(bkbn_atoms, i_gene->lookupRot->rotamerlet[h].coord, iRES->atom_ptr, i_gene->seq_position); energy_vdw = dummy_var_fix_vdw_energy.E_vdw; if(energy_vdw - var_var_ENERGy.E_vdw <= EPS) { best_rotamerlet = h; var_var_ENERGy.E_vdw = energy_vdw; } } if(iRES->num_methyl_groups != 0) { dummy_var_fix_vdw_energy = var_fixed_vdw_screen(bkbn_atoms, iROT->rotamerlet[best_rotamerlet].methyl_group->non_MeH, iRES->non_MeH_atom_ptr, i_gene->seq_position); var_var_ENERGy.E_vdw = dummy_var_fix_vdw_energy.E_vdw; for(x=1;x<=iRES->num_methyl_groups;++x) { dummy_vdw_energy1 = DBL_MAX; for(y=1;y<=3;++y) { dummy_var_fix_vdw_energy = var_fixed_vdw_screen(bkbn_atoms, iROT->rotamerlet[best_rotamerlet].methyl_group->MeH[x].coord[y], iRES->MeH_atom_ptr[x], i_gene->seq_position); energy_vdw = dummy_var_fix_vdw_energy.E_vdw; if(energy_vdw - dummy_vdw_energy1 <= EPS) { dummy_vdw_energy1 = energy_vdw; best_methyl = y; } } var_var_ENERGy.E_vdw += dummy_vdw_energy1; } } total_vdw_energy += var_var_ENERGy.E_vdw; } else { if(i_gene->choice_ptr->lookup_res_ptr->lookupRot[1].sideAtoms == NULL) { if(i_gene->choice_ptr->lookup_res_ptr->fixed_flag == 1) sprintf(structure_filename, "%s/coordinates/%d/%s.%d.structure.fixed", parameters->lookup_energy_table_directory, i_gene->seq_position, i_gene->choice_ptr->resparam_ptr->residuetype, i_gene->seq_position); else sprintf(structure_filename, "%s/coordinates/%d/%s.%d.structure", parameters->lookup_energy_table_directory, i_gene->seq_position, i_gene->choice_ptr->resparam_ptr->residuetype, i_gene->seq_position); load_coordinates_lookupRes(i_gene->choice_ptr->lookup_res_ptr, i_gene->choice_ptr->resparam_ptr, structure_filename); } q=1; while(i_gene->choice_ptr->resparam_ptr->atom_ptr[q]!=NULL) { sideAtoms[q].coord = i_gene->lookupRot->sideAtoms[q].coord; sideAtoms[q].born_radius = i_gene->lookupRot->sideAtoms[q].born_radius; sideAtoms[q].sasa = 0; sideAtoms[q].atom_ptr = i_gene->choice_ptr->resparam_ptr->atom_ptr[q]; sideAtoms[q].seq_position = i_gene->seq_position; ++q; } sideAtoms[q].seq_position = ENDFLAG; sideAtoms[q].sasa = ENDFLAG; sideAtoms[q].atom_ptr = NULL; var_fix_ENERGy = var_fixed_energy_calc(bkbn_atoms, sideAtoms,1); total_vdw_energy += var_fix_ENERGy.E_vdw; } j = i+1; chr->genes = chr->genes->nextgene; while(chr->genes->seq_position!=ENDFLAG) { if(chr->genes->choice_ptr->lookup_res_ptr->lookupRot[1].sideAtoms == NULL) { if(chr->genes->choice_ptr->lookup_res_ptr->fixed_flag == 1) sprintf(structure_filename, "%s/coordinates/%d/%s.%d.structure.fixed", parameters->lookup_energy_table_directory, chr->genes->seq_position, chr->genes->choice_ptr->resparam_ptr->residuetype, chr->genes->seq_position); else sprintf(structure_filename, "%s/coordinates/%d/%s.%d.structure", parameters->lookup_energy_table_directory, chr->genes->seq_position, chr->genes->choice_ptr->resparam_ptr->residuetype, chr->genes->seq_position); load_coordinates_lookupRes(chr->genes->choice_ptr->lookup_res_ptr, chr->genes->choice_ptr->resparam_ptr, structure_filename); } var_var_ENERGy = var_var_energy_calc(i_gene->lookupRot->sideAtoms, chr->genes->lookupRot->sideAtoms, i_gene->choice_ptr->resparam_ptr->atom_ptr, chr->genes->choice_ptr->resparam_ptr->atom_ptr, &do_these_hbond_flag, i_hbond_status, j_hbond_status); if(LOCAL_MINIMIZATION_FLAG == 1) { if(var_var_ENERGy.E_vdw!=0) { var_var_ENERGy.E_vdw = DBL_MAX; for(h=1;h<=iRES->rotamerlib_ptr->num_rotamerlets;++h) { for(k=1;k<=CHRRES->rotamerlib_ptr->num_rotamerlets;++k) { vdw_energy = var_var_vdw_screen_micro(i_gene->lookupRot->rotamerlet[h].coord, chr->genes->lookupRot->rotamerlet[k].coord, iRES->atom_ptr, CHRRES->atom_ptr); if(vdw_energy - var_var_ENERGy.E_vdw <= EPS) { best_rotamerlet_i = h; best_rotamerlet_j = k; var_var_ENERGy.E_vdw = vdw_energy; } } } /* methyl relaxation for best rotamerlet pair */ // if(MINIMIZE_METHYL_FLAG == 1) if(iRES->num_methyl_groups!=0 || CHRRES->num_methyl_groups!=0) /* relevant only if at least one of them has a methyl */ { if(iRES->num_methyl_groups!=0 && CHRRES->num_methyl_groups!=0) /* both have methyls */ { /* inter MeH energies */ for(x=1;x<=iRES->num_methyl_groups;++x) { for(y=1;y<=3;++y) { for(h=1;h<=CHRRES->num_methyl_groups;++h) { for(k=1;k<=3;++k) { MeH_MeH_energy[x][y][h][k] = var_var_vdw_screen_micro(i_gene->lookupRot->rotamerlet[best_rotamerlet_i].methyl_group->MeH[x].coord[y], chr->genes->lookupRot->rotamerlet[best_rotamerlet_j].methyl_group->MeH[h].coord[k], iRES->MeH_atom_ptr[x], CHRRES->MeH_atom_ptr[h]); } } } } /* MeH-non_MeH energies */ for(x=1;x<=iRES->num_methyl_groups;++x) { for(y=1;y<=3;++y) { MeH_nonMeH_energy_i[x][y] = var_var_vdw_screen_micro(i_gene->lookupRot->rotamerlet[best_rotamerlet_i].methyl_group->MeH[x].coord[y], chr->genes->lookupRot->rotamerlet[best_rotamerlet_j].methyl_group->non_MeH, iRES->MeH_atom_ptr[x], CHRRES->non_MeH_atom_ptr); } } for(h=1;h<=CHRRES->num_methyl_groups;++h) { for(k=1;k<=3;++k) { MeH_nonMeH_energy_j[h][k] = var_var_vdw_screen_micro(i_gene->lookupRot->rotamerlet[best_rotamerlet_i].methyl_group->non_MeH, chr->genes->lookupRot->rotamerlet[best_rotamerlet_j].methyl_group->MeH[h].coord[k], iRES->non_MeH_atom_ptr, CHRRES->MeH_atom_ptr[h]); } } /* nonMeH-non_MeH energies */ nonMeH_nonMeH_energy = var_var_vdw_screen_micro(i_gene->lookupRot->rotamerlet[best_rotamerlet_i].methyl_group->non_MeH, chr->genes->lookupRot->rotamerlet[best_rotamerlet_j].methyl_group->non_MeH, iRES->non_MeH_atom_ptr, CHRRES->non_MeH_atom_ptr); /* add up and find the best one */ var_var_ENERGy.E_vdw = DBL_MAX; for(conf_i=1; conf_i<=iRES->num_methyl_confs; ++conf_i) { for(conf_j=1; conf_j<=CHRRES->num_methyl_confs; ++conf_j) { vdw_energy = 0; for(x=1;x<=iRES->num_methyl_groups;++x) { for(h=1;h<=CHRRES->num_methyl_groups;++h) { vdw_energy += MeH_MeH_energy[x][BASE_THREE_CONVERT[conf_i][x]][h][BASE_THREE_CONVERT[conf_j][h]]; } vdw_energy += MeH_nonMeH_energy_i[x][BASE_THREE_CONVERT[conf_i][x]]; } for(h=1;h<=CHRRES->num_methyl_groups;++h) { vdw_energy += MeH_nonMeH_energy_j[h][BASE_THREE_CONVERT[conf_j][h]]; } if(vdw_energy - var_var_ENERGy.E_vdw <= EPS) { var_var_ENERGy.E_vdw = vdw_energy; } } } var_var_ENERGy.E_vdw += nonMeH_nonMeH_energy; } else /* only one has a methyl */ { if(iRES->num_methyl_groups!=0) { rotamerlet_Me = &i_gene->lookupRot->rotamerlet[best_rotamerlet_i]; resparam_ptr_Me = iRES; rotamerlet_non_Me = &chr->genes->lookupRot->rotamerlet[best_rotamerlet_j]; resparam_ptr_non_Me = CHRRES; } else { rotamerlet_Me = &chr->genes->lookupRot->rotamerlet[best_rotamerlet_j]; resparam_ptr_Me = CHRRES; rotamerlet_non_Me = &i_gene->lookupRot->rotamerlet[best_rotamerlet_i]; resparam_ptr_non_Me = iRES; } /* nonMeH-other_sidechain energies */ vdw_energy = var_var_vdw_screen_micro(rotamerlet_non_Me->coord, rotamerlet_Me->methyl_group->non_MeH, resparam_ptr_non_Me->atom_ptr, resparam_ptr_Me->non_MeH_atom_ptr); var_var_ENERGy.E_vdw = vdw_energy; /* MeH-other_sidechain energies */ for(h=1;h<=resparam_ptr_Me->num_methyl_groups;++h) { dummy_vdw_energy1 = DBL_MAX; for(k=1;k<=3;++k) { vdw_energy = var_var_vdw_screen_micro(rotamerlet_non_Me->coord, rotamerlet_Me->methyl_group->MeH[h].coord[k], resparam_ptr_non_Me->atom_ptr, resparam_ptr_Me->MeH_atom_ptr[h]); if(vdw_energy - dummy_vdw_energy1 <= EPS) { dummy_vdw_energy1 = vdw_energy; } } var_var_ENERGy.E_vdw += dummy_vdw_energy1; } } } } } total_vdw_energy += var_var_ENERGy.E_vdw; ++j; chr->genes = chr->genes->nextgene; } ++i; chr->genes = i_gene->nextgene; } chr->genes = chr->firstgene; free_memory(bkbn_atoms); free_memory(sideAtoms); free_memory(structure_filename); free_memory(i_hbond_status); free_memory(j_hbond_status); return(total_vdw_energy); }