/* EGAD: lookup_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 creating and de-allocating pair energy lookup tables */ #include "lookup_table.h" float NON_INTERACT_PTR=0; float FIXED_POSITION_PTR=0; LOOKUP_ENERGY_ROTAMER_X *NON_INTERACT_LOOKUP_ROT_X=NULL; LOOKUP_ENERGY_RESIDUE_X *NON_INTERACT_LOOKUP_RES_X=NULL; int **BASE_THREE_CONVERT=NULL; bool *FIXED_HBOND_SATISF_ARRAY; ATOMRESPARAM **FIXED_ATOM_HBOND_ATOMRESPARAM; int FIXED_ATOM_HBOND_STATUS_ARRAYLENGTH, BKBN_HBOND_STATUS_ARRAYLENGTH; /* initialize non-interaction pointers for the lookuptable, and BASE_THREE_CONVERT matrix for local minimization of methyl groups. called by generate_lookup_table */ void initialize_lookuptable_pointers() { int i, j; double base3_exp, base3_diff, base3_quot; extern float NON_INTERACT_PTR, FIXED_POSITION_PTR; extern LOOKUP_ENERGY_ROTAMER_X *NON_INTERACT_LOOKUP_ROT_X; extern LOOKUP_ENERGY_RESIDUE_X *NON_INTERACT_LOOKUP_RES_X; extern int **BASE_THREE_CONVERT; /* some dummy pointers for the lookup table */ /* not really pointers, but these values are unlikely to ever come up in an actual run */ NON_INTERACT_PTR = -1234.5; FIXED_POSITION_PTR = -2345.6; NON_INTERACT_LOOKUP_RES_X = (LOOKUP_ENERGY_RESIDUE_X *)malloc(sizeof(LOOKUP_ENERGY_RESIDUE_X)); NON_INTERACT_LOOKUP_ROT_X = (LOOKUP_ENERGY_ROTAMER_X *)malloc(sizeof(LOOKUP_ENERGY_ROTAMER_X)); /* used for generating conformations for local minimization of methyl groups */ BASE_THREE_CONVERT = (int **)calloc((int)(pow(3, MAX_NUM_METHYL) + 1), sizeof(int *)); for(i=1;i<=pow(3, MAX_NUM_METHYL);++i) BASE_THREE_CONVERT[i] = (int *)calloc(MAX_NUM_METHYL + 1, sizeof(int)); for(i=1;i<=pow(3, MAX_NUM_METHYL);++i) { base3_diff = i-1.0; for(j=MAX_NUM_METHYL;j>=1;--j) { base3_exp = (double)pow(3,(j-1)); base3_quot = base3_diff/base3_exp; if(base3_quot<1) BASE_THREE_CONVERT[i][j] = 1; else { if(base3_quot<2) BASE_THREE_CONVERT[i][j] = 2; else BASE_THREE_CONVERT[i][j] = 3; base3_diff = base3_diff - (BASE_THREE_CONVERT[i][j]-1)*base3_exp; } } } } void free_LOOKUP_ROTAMER_X_HBOND(LOOKUP_ROTAMER_X_HBOND **lookupRotX_hbond) { int i; if((*lookupRotX_hbond) == NULL) return; if((*lookupRotX_hbond)->hbond_rot_pair == NULL) return; for(i=0;i<(*lookupRotX_hbond)->num_hbond_rot_pairs;++i) { free_memory((*lookupRotX_hbond)->hbond_rot_pair[i].i_hbond_status); free_memory((*lookupRotX_hbond)->hbond_rot_pair[i].j_hbond_status); } free_memory((*lookupRotX_hbond)->hbond_rot_pair); free_memory((*lookupRotX_hbond)); } /* This function de-allocates memory occupied by the lookup table. Called by input_stuff.cpp: free_PROTEIN */ void free_lookup_table(LOOKUP_ENERGY *lookupEnergy, VARIABLE_POSITION *varPos) { int i, j; /* variable position counters */ int i_res, j_res; /* variable residuetype counters */ int i_res_rot; /* variable rotamer counters */ int numVarPositions; extern LOOKUP_ENERGY_ROTAMER_X *NON_INTERACT_LOOKUP_ROT_X; extern LOOKUP_ENERGY_RESIDUE_X *NON_INTERACT_LOOKUP_RES_X; if(lookupEnergy==NULL) return; i=1; while(varPos[i].seq_position!=ENDFLAG) ++i; numVarPositions = i-1; for(i=1;i<=numVarPositions;++i) { for(i_res=1; i_res<=varPos[i].number_of_choices; ++i_res) { if(lookupEnergy[i].lookupRes[i_res].lookupResRes!=NULL) { free_memory(lookupEnergy[i].lookupRes[i_res].lookupResRes); lookupEnergy[i].lookupRes[i_res].lookupResRes = NULL; } for(i_res_rot=1; i_res_rot<= lookupEnergy[i].lookupRes[i_res].sizeof_lookupRot; ++i_res_rot) { if(ERESROT.P != NULL) free_memory(ERESROT.P); if(ERESROT.sidechain_hbond_status != NULL) free_memory(ERESROT.sidechain_hbond_status); if(ERESROT.bkbn_hbond_indicies != NULL) free_memory(ERESROT.bkbn_hbond_indicies); if(ERESROT.lookupX!=NULL) { for(j=(i+1);j<=numVarPositions;++j) { if(ERESROT.lookupX[j-i].lookupResX!=NULL && ERESROT.lookupX[j-i].lookupResX!=NON_INTERACT_LOOKUP_RES_X) { for(j_res=0;j_reschoice[i_res].resparam_ptr->num_hbonding_sidechain_atoms!=0) { if(varpos->core_flag == 's') return(varpos->choice[i_res].resparam_ptr->E_specificity_surface); return(0.0); } if(varpos->core_flag == 's') return(varpos->choice[i_res].resparam_ptr->E_specificity_surface); if(varpos->core_flag == 'i') return(varpos->choice[i_res].resparam_ptr->E_specificity_interfacial); return(varpos->choice[i_res].resparam_ptr->E_specificity_core); // core } /* structure used for a temporary holding variable in generate_lookup_table */ typedef struct { ROTAMER *rotamer; /* pointer to the rotamer */ int keep_flag; /* 1 = keep this rotamer */ int index; /* rotamer index */ var_fix_ENERGY energy_var_fix; double E_forcefield; double E_gbsa; double P; /* boltzmann probab of this rotamer based on energy */ double sasa_total; double E_asp; double sasa_hydrophobic; mini_pdbATOM *sideAtoms; /* mini_pdbATOM array rotamer p at this position i */ /* SASA values are the areas exposed on this sidechain in the presence of fixedatoms and no other sidechains */ /* sidechain Born radii listed for charged atoms */ ROTAMERLET *rotamerlet; /* rotamerlet conformations for local minimization */ bool *bkbn_hbond_status; bool *sidechain_hbond_status; } PSEUDO_LOOKUP_ENERGY_ROTAMER; /* This function generates a lookup table for all possible var_residue-fixed_atom and var-var pairwise interaction energies for fixed backbone optimizations. Links protein->var_pos lookup_table pointers to appropriate locations. Assumes that protein has been initialized with input_stuff.cpp: input_stuff All the rotamer optimization jobs assume that this function has been called. */ void generate_lookup_table(PROTEIN *protein) { int i, j; /* variable position counters */ int i_res, j_res; /* variable residuetype counters */ int i_res_rot, j_res_rot, numRot, numRot0; /* variable rotamer counters */ CHROMOSOME chr_i, *chr; int N_position, N_position_pseudo, N_position_bkbn_hbond_status; int dummy; int another_flag, proline_flag, end_of_sidechain_index, go_flag; double alpha, f; var_fix_vdw_ENERGY energy_rotamerlet, dummy_energy_rotamerlet; var_fix_vdw_ENERGY dummy_energy_Me, energy_Me; double E_sasa_pro_gly, sasa_pro_gly_total, one_over_numRot, E_sasa_pro_gly_asp, sasa_pro_gly_hydrophobic; double denom, E_temp; double E_pol_intra; int n, q, h, k, z; mini_pdbATOM *dummyminipdb = NULL, *side_minipdb = NULL; mini_pdbATOM *fixed_atoms = NULL; mini_pdbATOM *actual_CB_atom = NULL; mini_pdbATOM *pseudo_sidechain = NULL; mini_pdbATOM *pseudo_fixed_atoms = NULL, *holding_pseudo_fixed_atoms = NULL; int *pro_and_gly_pos = NULL, num_pro_gly; BACKBONE *bkbn = NULL; SIDECHAIN side; int num_heavy_fixedplus1; PSEUDO_LOOKUP_ENERGY_ROTAMER *pseudo_lookupRotamer = NULL; ROTAMERLIB *rotamerlib_ptr = NULL; int best_methyl; int best_rotamer, best_rotamerlet; int x, y, numVarPositions, calc_var_fix_energy_flag, calc_var_var_flag; CARTESIAN **foundation_coords = NULL; double dihedangle; char *dummystring, *directoryname; milli_pdbATOM **i_atoms = NULL, **j_atoms = NULL; double xmax, ymax, zmax, xmin, ymin, zmin; double E_forcefield_min, E_gbsa_min; double num_rot_position, EE; int fixed_positions_flag; int ii, jj, end_real_fixed_atoms_index; int last_seq_position, first_position; short int *holding_shortint_array; time_t now; double boltz_denom; ATOMRESPARAM **holding_atomresparam; extern double **CHARGE_PRODUCT; extern float NON_INTERACT_PTR; extern LOOKUP_ENERGY_ROTAMER_X *NON_INTERACT_LOOKUP_ROT_X; extern LOOKUP_ENERGY_RESIDUE_X *NON_INTERACT_LOOKUP_RES_X; extern ATOMRESPARAM PSEUDO_SIDECHAIN_ATOMRESPARAM, BKBN_N_ATOMRESPARAM,BKBN_CA_ATOMRESPARAM, BKBN_C_ATOMRESPARAM,BKBN_O_ATOMRESPARAM, BKBN_HA_ATOMRESPARAM; extern CARTESIAN_VECTOR ***VECTOR_PTR; extern int MAX_SEQ_POSITION; double rotamer_complexity; extern int MAX_ROTAMERS, MAX_MUTATIONS; extern int LOCAL_MINIMIZATION_FLAG,LOGFILE_FLAG; extern int SASA_FLAG, GB_FLAG, COULOMB_FLAG, PAIR_ENERGY_TABLE_FLAG, TORSION_FLAG; extern int GET_PID; extern char *CURRENT_WORKING_DIRECTORY, *LOOKUP_TABLE_DIRECTORY; LOOKUP_ENERGY *lookupEnergy; VARIABLE_POSITION *varPos; PARAMETERS *parameters; mini_pdbATOM *input_fixed_atoms; FILE *logfile_ptr; MENDEL i_gene, j_gene; bool *fixed_atoms_hbond_satisfaction_array, *holding_fixed_hbond_satisf_array, free_ligand_flag; extern bool *FIXED_HBOND_SATISF_ARRAY; int num_fixed_atom_hbonding_atoms, new_num_fixed_atom_hbonding_atoms; LOOKUP_ROTAMER_X_HBOND *lookupRotX_hbond; double dihed_CB; /* for building pse along CB */ /* allocate memory for SIDECHAIN side.atom array */ side.atom=(micro_pdbATOM *)calloc(MAX_RES_SIZE,sizeof(micro_pdbATOM)); strcpy(LOOKUP_TABLE_DIRECTORY,protein->parameters.lookup_energy_table_directory); logfile_ptr=NULL; input_fixed_atoms = protein->mini_fixed_atoms; varPos = protein->var_pos; parameters = &protein->parameters; protein->lookupEnergy = (LOOKUP_ENERGY *)calloc((parameters->numVarPositions+2),sizeof(LOOKUP_ENERGY)); lookupEnergy = protein->lookupEnergy; /* if the lookuptable exists on disk, make sure that the forcefield and backbone are identical; if it doesn't exit, create the directory and a file with the forcefield and template info */ check_lookup_table(protein); j=1;j_res=1;j_res_rot=1; best_rotamer = ENDFLAG; best_rotamerlet = ENDFLAG; /* Find out if any of the positions are defined by the user as fixed; if so, fixed_positions_flag = 1 and lookupResRes containing interaction lists are NOT saved to the lookuptable on disk for these */ dummystring = (char *)calloc(MAXLINE, sizeof(char)); directoryname = (char *)calloc(MXLINE_INPUT, sizeof(char)); fixed_positions_flag = 0; i=1; while(varPos[i].seq_position!=ENDFLAG) { if(varPos[i].fixed_flag == 1) fixed_positions_flag = 1; ++i; } numVarPositions = i-1; ii=1; num_fixed_atom_hbonding_atoms=0; first_position = input_fixed_atoms[ii].seq_position; while(input_fixed_atoms[ii].seq_position!=ENDFLAG) { last_seq_position = input_fixed_atoms[ii].seq_position; if(input_fixed_atoms[ii].atom_ptr->donorflag[0] !=' ') ++num_fixed_atom_hbonding_atoms; ++ii; } fixed_atoms_hbond_satisfaction_array = (bool *)calloc(num_fixed_atom_hbonding_atoms+2,sizeof(bool)); falsify_bool_array(fixed_atoms_hbond_satisfaction_array,1,num_fixed_atom_hbonding_atoms); FIXED_ATOM_HBOND_ATOMRESPARAM = (ATOMRESPARAM **)calloc(num_fixed_atom_hbonding_atoms+2,sizeof(ATOMRESPARAM *)); ii=1; x=0; while(input_fixed_atoms[ii].seq_position!=ENDFLAG) { if(input_fixed_atoms[ii].atom_ptr->donorflag[0] !=' ') { ++x; fixed_atoms_hbond_satisfaction_array[x] = input_fixed_atoms[ii].hbond_satisfied; FIXED_ATOM_HBOND_ATOMRESPARAM[x] = input_fixed_atoms[ii].atom_ptr; } ++ii; } /* printf("before pro/gly: %d ",count_false(fixed_atoms_hbond_satisfaction_array,1,num_fixed_atom_hbonding_atoms)); fprintf_bool_array(stdout,fixed_atoms_hbond_satisfaction_array,1,num_fixed_atom_hbonding_atoms); */ parameters->numVarPositions = numVarPositions; /* intialize these variables only the first time this function is called */ if(NON_INTERACT_PTR==0) initialize_lookuptable_pointers(); calc_var_fix_energy_flag = 1; /* create the "var_fix" directory in the lookup_table directory path if it doesn't exist */ sprintf(directoryname,"var_fix"); while(create_directory(directoryname, parameters->lookup_energy_table_directory)==0); /* load lookupRes; ie: var-fix energies; returns 0 if all the energies were loaded, returns 1 if one or more energies need to be calculated */ calc_var_fix_energy_flag = load_lookupRes_from_disk(lookupEnergy, varPos, num_fixed_atom_hbonding_atoms, parameters); if(calc_var_fix_energy_flag == 1) /* need to create var-fix energies, etc for at least one residue in lookup table */ { /* create the "coordinates" directory in the lookup_table directory path if it doesn't exist */ sprintf(directoryname,"coordinates"); while(create_directory(directoryname, parameters->lookup_energy_table_directory)==0); /* allocate some memory */ pseudo_lookupRotamer = (PSEUDO_LOOKUP_ENERGY_ROTAMER *)calloc(MAX_ROTAMERS, sizeof(PSEUDO_LOOKUP_ENERGY_ROTAMER)); bkbn = (BACKBONE *)calloc(MAX_RESIDUES, sizeof(BACKBONE)); actual_CB_atom = (mini_pdbATOM *)calloc(MAX_SEQ_POSITION + 10, sizeof(mini_pdbATOM)); foundation_coords = (CARTESIAN **)calloc(4, sizeof(CARTESIAN *)); pro_and_gly_pos = (int *)calloc(MAX_RESIDUES, sizeof(int)); dummyminipdb = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM)); side_minipdb = (mini_pdbATOM *)calloc(MAX_RES_SIZE, sizeof(mini_pdbATOM)); holding_pseudo_fixed_atoms = (mini_pdbATOM *)calloc(MAX_RES_SIZE, sizeof(mini_pdbATOM)); /* create chr_i CHROMOSOME */ chr_i.genes = (MENDEL *)malloc(sizeof(MENDEL)); /* gene i */ chr_i.firstgene = chr_i.genes; chr_i.genes->nextgene = (MENDEL *)malloc(sizeof(MENDEL)); chr_i.genes = chr_i.genes->nextgene; /* end gene */ chr_i.genes->seq_position = ENDFLAG; chr_i.genes=chr_i.firstgene; fixed_atoms = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM)); fixed_atoms[0].sasa=0; fixed_atoms[0].seq_position=0; actual_CB_atom[0].sasa=0; actual_CB_atom[0].seq_position=0; pseudo_sidechain = (mini_pdbATOM *)calloc(MAX_RESIDUES, sizeof(mini_pdbATOM)); pseudo_sidechain[0].sasa=0; pseudo_sidechain[0].seq_position=0; pseudo_fixed_atoms = (mini_pdbATOM *)calloc((MAX_ATOMS), sizeof(mini_pdbATOM)); pseudo_fixed_atoms[0].sasa=0; pseudo_fixed_atoms[0].seq_position=0; extract_mini_to_bkbn(input_fixed_atoms, bkbn); i=1; j=1; k=1; n=1; q=1; while(input_fixed_atoms[i].seq_position!=ENDFLAG) /* fixed_atoms in this function doesn't include CB */ { if(strcmp(input_fixed_atoms[i].atom_ptr->atomname, "CB")==0 || strcmp(input_fixed_atoms[i].atom_ptr->atomname, "2HA")==0) { actual_CB_atom[input_fixed_atoms[i].seq_position] = input_fixed_atoms[i]; /* create pseudoatoms */ /* original pse build - builds using this torsional description H->N->CA->PSE; as pointed out by Mark Voorhies, the location of PSE w/ respect to CA-CB vector will change w/ phi/psi */ /* pseudo_sidechain[k].coord = chi2xyz(&bkbn[k].CA, &bkbn[k].N, &bkbn[k].H, &PSEUDO_SIDECHAIN_ATOMRESPARAM.bondlength, &PSEUDO_SIDECHAIN_ATOMRESPARAM.bondangle, &PSEUDO_SIDECHAIN_ATOMRESPARAM.dihedral); */ /* build pse along CB */ dihed_CB = bkbn[k].phi-120; pseudo_sidechain[k].coord = chi2xyz(&bkbn[k].CA, &bkbn[k].N, &bkbn[k-1].C, &PSEUDO_SIDECHAIN_ATOMRESPARAM.bondlength, &PSEUDO_SIDECHAIN_ATOMRESPARAM.bondangle, &dihed_CB); pseudo_sidechain[k].seq_position = actual_CB_atom[input_fixed_atoms[i].seq_position].seq_position; pseudo_sidechain[k].box=0; pseudo_sidechain[k].atom_ptr = &PSEUDO_SIDECHAIN_ATOMRESPARAM; // U atom for CB; for ligands, the backbone doesn't really exist, so don't include a pseudosidechain if(actual_CB_atom[input_fixed_atoms[i].seq_position].atom_ptr->atom_ptr->atom_class == ENDFLAG) pseudo_sidechain[k].atom_ptr = actual_CB_atom[input_fixed_atoms[i].seq_position].atom_ptr; ++k; } else { fixed_atoms[j] = input_fixed_atoms[i]; // if(fixed_atoms[j].atom_ptr->atom_ptr->sasa_radius>0) { pseudo_fixed_atoms[n] = fixed_atoms[j]; pseudo_fixed_atoms[n].atom_ptr = (ATOMRESPARAM *)malloc(sizeof(ATOMRESPARAM)); *pseudo_fixed_atoms[n].atom_ptr = *fixed_atoms[j].atom_ptr; pseudo_fixed_atoms[n].atom_ptr->atom_ptr = (ATOMPARAM *)malloc(sizeof(ATOMPARAM)); *pseudo_fixed_atoms[n].atom_ptr->atom_ptr = *fixed_atoms[j].atom_ptr->atom_ptr; if(strcmp(pseudo_fixed_atoms[n].atom_ptr->atomname, "N")==0) { *pseudo_fixed_atoms[n].atom_ptr->atom_ptr = *BKBN_N_ATOMRESPARAM.atom_ptr; } else if(strcmp(pseudo_fixed_atoms[n].atom_ptr->atomname, "C")==0) { *pseudo_fixed_atoms[n].atom_ptr->atom_ptr = *BKBN_C_ATOMRESPARAM.atom_ptr; } else if(strcmp(pseudo_fixed_atoms[n].atom_ptr->atomname, "CA")==0) { *pseudo_fixed_atoms[n].atom_ptr->atom_ptr = *BKBN_CA_ATOMRESPARAM.atom_ptr; } else if(strcmp(pseudo_fixed_atoms[n].atom_ptr->atomname, "O")==0) { *pseudo_fixed_atoms[n].atom_ptr->atom_ptr = *BKBN_O_ATOMRESPARAM.atom_ptr; } else if(strcmp(pseudo_fixed_atoms[n].atom_ptr->atomname, "HA")==0) { *pseudo_fixed_atoms[n].atom_ptr->atom_ptr = *BKBN_HA_ATOMRESPARAM.atom_ptr; } ++n; } ++j; } ++i; } end_real_fixed_atoms_index = j; /* actual_CB_atom[k].seq_position=ENDFLAG; actual_CB_atom[k].sasa = ENDFLAG; */ pseudo_sidechain[k].seq_position=ENDFLAG; pseudo_sidechain[k].sasa = ENDFLAG; pseudo_sidechain[k].atom_ptr = NULL; pseudo_fixed_atoms[n].seq_position=ENDFLAG; pseudo_fixed_atoms[n].atom_ptr=NULL; pseudo_fixed_atoms[n].sasa = ENDFLAG; k=1; num_pro_gly=0; while(pseudo_sidechain[k].seq_position!=ENDFLAG) { i=1; another_flag=0; while(i<=numVarPositions && another_flag==0) /* find variable positions */ { if(pseudo_sidechain[k].seq_position == varPos[i].seq_position) another_flag=1; ++i; } if(another_flag==0) /* not a variable position; is a pro or gly */ { ++num_pro_gly; pro_and_gly_pos[num_pro_gly] = pseudo_sidechain[k].seq_position; } ++k; } /* calculate SASA for pro/gly positions */ N_position_pseudo=1; N_position=1; E_sasa_pro_gly = 0; E_sasa_pro_gly_asp=0; N_position_bkbn_hbond_status=1; sasa_pro_gly_total=0; sasa_pro_gly_hydrophobic = 0; for(i=1;i<=num_pro_gly;++i) { fixed_atoms[j] = actual_CB_atom[pro_and_gly_pos[i]]; ++j; while(pseudo_fixed_atoms[N_position_pseudo].seq_position!=pro_and_gly_pos[i]) ++N_position_pseudo; while(fixed_atoms[N_position].seq_position!=pro_and_gly_pos[i]) { if(fixed_atoms[N_position].atom_ptr->donorflag[0] != ' ') ++N_position_bkbn_hbond_status; ++N_position; } proline_flag=0; q=N_position; n=1; h=N_position_pseudo; /* for position i, replace the pseudobackbone atoms w/ the real thing */ while(fixed_atoms[q].seq_position==pro_and_gly_pos[i]) /* store the pseudoatoms for pos i in holding_pseudo_fixed_atoms */ { if( /* fixed_atoms[q].atom_ptr->atom_ptr->sasa_radius>0 && */ fixed_atoms[q].atom_ptr->atomnumber != 8) { holding_pseudo_fixed_atoms[n] = pseudo_fixed_atoms[h]; pseudo_fixed_atoms[h] = fixed_atoms[q]; ++n; ++h; } if(strcmp(fixed_atoms[q].atom_ptr->atomname,"CD")==0) /* is a proline */ proline_flag=1; ++q; } holding_pseudo_fixed_atoms[n].sasa=ENDFLAG; holding_pseudo_fixed_atoms[n].seq_position = ENDFLAG; holding_pseudo_fixed_atoms[n].atom_ptr = NULL; dummyminipdb[0].seq_position=ENDFLAG; dummyminipdb[0].atom_ptr=NULL; q=1; while(pseudo_fixed_atoms[q].seq_position!=ENDFLAG) { dummyminipdb[q] = pseudo_fixed_atoms[q]; ++q; } h=1; while(pseudo_sidechain[h].seq_position!=ENDFLAG) { if(pseudo_sidechain[h].seq_position!=pro_and_gly_pos[i]) { dummyminipdb[q] = pseudo_sidechain[h]; ++q; } ++h; } // if(actual_CB_atom[pro_and_gly_pos[i]].atom_ptr->atom_ptr->sasa_radius>0) { dummyminipdb[q] = actual_CB_atom[pro_and_gly_pos[i]]; ++q; } dummyminipdb[q].sasa=ENDFLAG; dummyminipdb[q].seq_position=ENDFLAG; if(SASA_FLAG!=0) { h=N_position_pseudo; while(dummyminipdb[h].seq_position == dummyminipdb[N_position_pseudo].seq_position) { dummyminipdb[h].hbond_satisfied=0; ++h; } SASAcalc(dummyminipdb); exposed_hbond_groups(dummyminipdb); h=N_position_pseudo; jj=N_position_bkbn_hbond_status; while(dummyminipdb[h].seq_position == dummyminipdb[N_position_pseudo].seq_position) { sasa_pro_gly_total += dummyminipdb[h].sasa; E_sasa_pro_gly_asp += dummyminipdb[h].sasa*dummyminipdb[h].atom_ptr->atom_ptr->asp; if(dummyminipdb[h].atom_ptr->atom_ptr->atomtype[0] == 'C' || dummyminipdb[h].atom_ptr->atom_ptr->atomtype[0] == 'S' || dummyminipdb[h].atom_ptr->atom_ptr->atomtype[0] == 'H') sasa_pro_gly_hydrophobic += dummyminipdb[h].sasa; if(dummyminipdb[h].atom_ptr->donorflag[0] != ' ') { fixed_atoms_hbond_satisfaction_array[jj] = dummyminipdb[h].hbond_satisfied || fixed_atoms_hbond_satisfaction_array[jj]; ++jj; } ++h; } /* CB sasa */ h=q-1; sasa_pro_gly_total += dummyminipdb[h].sasa; E_sasa_pro_gly_asp += dummyminipdb[h].sasa*dummyminipdb[h].atom_ptr->atom_ptr->asp; if(dummyminipdb[h].atom_ptr->atom_ptr->atomtype[0] == 'C' || dummyminipdb[h].atom_ptr->atom_ptr->atomtype[0] == 'H' || dummyminipdb[h].atom_ptr->atom_ptr->atomtype[0] == 'S') sasa_pro_gly_hydrophobic += dummyminipdb[h].sasa; E_sasa_pro_gly += GENERAL_ASP*sasa_pro_gly_total; } n=1; q=N_position_pseudo; /* replace backbone atoms w/ pseudoatoms in pseudo_fixed_atoms */ while(holding_pseudo_fixed_atoms[n].seq_position!=ENDFLAG) { pseudo_fixed_atoms[q] = holding_pseudo_fixed_atoms[n]; ++n; ++q; } } q=1; while(pseudo_sidechain[q].seq_position != ENDFLAG) { fixed_atoms[j] = pseudo_sidechain[q]; ++q; ++j; } fixed_atoms[j].seq_position=ENDFLAG; fixed_atoms[j].sasa = ENDFLAG; N_position_pseudo=1; N_position=1; N_position_bkbn_hbond_status=1; VECTOR_PTR = NULL; /* printf("after pro/gly: %d ",count_false(fixed_atoms_hbond_satisfaction_array,1,num_fixed_atom_hbonding_atoms)); fprintf_bool_array(stdout,fixed_atoms_hbond_satisfaction_array,1,num_fixed_atom_hbonding_atoms); */ /* start actual var-fix pair energy calculations */ for(i=1;i<=numVarPositions; ++i) /* position i loop */ { chr_i.genes = chr_i.firstgene; if(varPos[i].lookup_energy_ptr == NULL) { varPos[i].lookup_energy_ptr = &lookupEnergy[i]; /* copy lookup_energy_ptr */ lookupEnergy[i].seq_position = varPos[i].seq_position; /* copying sequence position */ /* allocate memory for LOOKUP_ENERGY_RESIDUE */ lookupEnergy[i].sizeof_lookupRes = varPos[i].number_of_choices; lookupEnergy[i].lookupRes = (LOOKUP_ENERGY_RESIDUE *)calloc((lookupEnergy[i].sizeof_lookupRes+1),sizeof(LOOKUP_ENERGY_RESIDUE)); } chr_i.genes->seq_position = varPos[i].seq_position; chr_i.genes->varpos_ptr = &varPos[i]; while(pseudo_fixed_atoms[N_position_pseudo].seq_position!=varPos[i].seq_position) { ++N_position_pseudo; } while(fixed_atoms[N_position].seq_position!=varPos[i].seq_position) { if(fixed_atoms[N_position].atom_ptr->donorflag[0] != ' ') ++N_position_bkbn_hbond_status; ++N_position; } /* generate dummyminipdb array which is used for SASA calcs; includes pseudo backbone and sidechain atoms for all positions except i; real backbone atoms for i. identify index N_position_pseudo of first backbone atom for this residue in this array */ q=N_position; n=1; h=N_position_pseudo; /* for position i, replace the pseudobackbone atoms w/ the real thing */ while(fixed_atoms[q].seq_position==varPos[i].seq_position) /* store the pseudoatoms for pos i in holding_pseudo_fixed_atoms */ { // if(fixed_atoms[q].atom_ptr->atom_ptr->sasa_radius>0) { holding_pseudo_fixed_atoms[n] = pseudo_fixed_atoms[h]; pseudo_fixed_atoms[h] = fixed_atoms[q]; ++n; ++h; } ++q; } holding_pseudo_fixed_atoms[n].sasa=ENDFLAG; holding_pseudo_fixed_atoms[n].seq_position = ENDFLAG; q=1; while(pseudo_fixed_atoms[q].seq_position!=ENDFLAG) { dummyminipdb[q] = pseudo_fixed_atoms[q]; ++q; } h=1; while(pseudo_sidechain[h].seq_position!=ENDFLAG) { if(pseudo_sidechain[h].seq_position!=varPos[i].seq_position) { dummyminipdb[q] = pseudo_sidechain[h]; ++q; } ++h; } num_heavy_fixedplus1 = q; for(i_res=1; i_res<=varPos[i].number_of_choices; ++i_res) /* residuetype i_res loop */ { if(varPos[i].choice[i_res].lookup_res_ptr == NULL) /* hasn't been calculated */ { varPos[i].choice[i_res].bkbn = &varPos[i].bkbn; varPos[i].choice[i_res].lookup_res_ptr = &lookupEnergy[i].lookupRes[i_res]; chr_i.genes->choice_ptr = &varPos[i].choice[i_res]; chr_i.genes->j_choice_index = i_res-1; lookupEnergy[i].lookupRes[i_res].residuecode = varPos[i].choice[i_res].resparam_ptr->residuecode; actual_CB_atom[varPos[i].seq_position].atom_ptr = &varPos[i].choice[i_res].resparam_ptr->atom[8]; /* assign proper ptr to CB */ lookupEnergy[i].lookupRes[i_res].lookupResRes = NULL; /* attach correct atom_ptr for backbone atoms for this residue */ q=N_position_pseudo; /* N_position_pseudo = index of first bkbn atom for this residue in dummyminipdb */ while(dummyminipdb[q].seq_position == dummyminipdb[N_position_pseudo].seq_position) { dummyminipdb[q].atom_ptr = &varPos[i].choice[i_res].resparam_ptr->sorted_atom[search_atomname_ATOMRESPARAM(dummyminipdb[q].atom_ptr->atomname, varPos[i].choice[i_res].resparam_ptr->sorted_atom)]; ++q; } numRot = 0; E_forcefield_min = DBL_MAX; E_gbsa_min=DBL_MAX; /* rotamer i_res_rot loop */ for(i_res_rot=1; i_res_rot<= varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot) { pseudo_lookupRotamer[i_res_rot].index = i_res_rot; pseudo_lookupRotamer[i_res_rot].bkbn_hbond_status = NULL; pseudo_lookupRotamer[i_res_rot].sidechain_hbond_status = NULL; /* copy chi angles */ chr_i.genes->chi = varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot].chi; chr_i.genes->lookupRot_index = i_res_rot; /* powell-minimize the rotamer (with each dihedral within 0.5 sigma of library values) */ if(LOCAL_MINIMIZATION_FLAG == 2) { minimize_rotamer_on_backbone(protein, chr_i.genes, fixed_atoms, &actual_CB_atom[varPos[i].seq_position]); } pseudo_lookupRotamer[i_res_rot].rotamer = (ROTAMER *)malloc(sizeof(ROTAMER)); *pseudo_lookupRotamer[i_res_rot].rotamer = varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot]; /* calculate sidechain structure */ pseudo_lookupRotamer[i_res_rot].sideAtoms = (mini_pdbATOM *)calloc(varPos[i].choice[i_res].resparam_ptr->numberofAtoms+5, sizeof(mini_pdbATOM)); GENE_to_mini_pdbATOM(&chr_i.genes,side_minipdb); pseudo_lookupRotamer[i_res_rot].sideAtoms[1] = actual_CB_atom[varPos[i].seq_position]; h=2; while(side_minipdb[h-1].seq_position!=ENDFLAG) { pseudo_lookupRotamer[i_res_rot].sideAtoms[h] = side_minipdb[h-1]; ++h; } pseudo_lookupRotamer[i_res_rot].sideAtoms[h].sasa=ENDFLAG; pseudo_lookupRotamer[i_res_rot].sideAtoms[h].seq_position=ENDFLAG; pseudo_lookupRotamer[i_res_rot].sideAtoms[h].atom_ptr=NULL; pseudo_lookupRotamer[i_res_rot].sideAtoms[0].sasa=ENDFLAG; pseudo_lookupRotamer[i_res_rot].sideAtoms[0].seq_position=ENDFLAG; pseudo_lookupRotamer[i_res_rot].sideAtoms[0].atom_ptr=NULL; /* append backbone atoms for i to the end of the pseudo_lookupRotamer[i_res_rot].sideAtoms array */ ++h; /* move past the ENDFLAG atom_ptr=NULL end of array marker */ q=N_position_pseudo; /* N_position_pseudo = index of first bkbn atom for this residue in dummyminipdb */ while(dummyminipdb[q].seq_position == dummyminipdb[N_position_pseudo].seq_position) { pseudo_lookupRotamer[i_res_rot].sideAtoms[h] = dummyminipdb[q]; ++h; ++q; } pseudo_lookupRotamer[i_res_rot].sideAtoms[h].seq_position=ENDFLAG2; pseudo_lookupRotamer[i_res_rot].sideAtoms[h].atom_ptr=NULL; /* put born_bonded into pseudo_lookupRotamer[i_res_rot].sideAtoms[].born_radius */ if(GB_FLAG==1) { q=1; while(pseudo_lookupRotamer[i_res_rot].sideAtoms[q].seq_position!=ENDFLAG) { if(pseudo_lookupRotamer[i_res_rot].sideAtoms[q].atom_ptr->charge!=0) if(pseudo_lookupRotamer[i_res_rot].sideAtoms[q].atom_ptr->atom_ptr->atom_class!=ENDFLAG) pseudo_lookupRotamer[i_res_rot].sideAtoms[q].born_radius = VRESROT.born_factor[q]; ++q; } } /* calculate energy_var_fix and born radii for this sidechain/rotamer */ pseudo_lookupRotamer[i_res_rot].energy_var_fix = var_fixed_energy_calc(fixed_atoms, pseudo_lookupRotamer[i_res_rot].sideAtoms,1); /* get hydrogen-bonding statuses */ pseudo_lookupRotamer[i_res_rot].bkbn_hbond_status = (bool *)calloc(num_fixed_atom_hbonding_atoms+2,sizeof(bool)); falsify_bool_array(pseudo_lookupRotamer[i_res_rot].bkbn_hbond_status,1,num_fixed_atom_hbonding_atoms); ii=1; jj=1; while(fixed_atoms[ii].seq_position!=ENDFLAG) { if(fixed_atoms[ii].atom_ptr->donorflag[0] != ' ') { pseudo_lookupRotamer[i_res_rot].bkbn_hbond_status[jj] = fixed_atoms[ii].hbond_satisfied || fixed_atoms_hbond_satisfaction_array[jj]; ++jj; } ++ii; } if(varPos[i].choice[i_res].resparam_ptr->num_hbonding_sidechain_atoms!=0) { pseudo_lookupRotamer[i_res_rot].sidechain_hbond_status = (bool *)calloc(varPos[i].choice[i_res].resparam_ptr->num_hbonding_sidechain_atoms+2,sizeof(bool)); falsify_bool_array(pseudo_lookupRotamer[i_res_rot].sidechain_hbond_status,1, varPos[i].choice[i_res].resparam_ptr->num_hbonding_sidechain_atoms); ii=1; jj=1; while(pseudo_lookupRotamer[i_res_rot].sideAtoms[ii].seq_position!=ENDFLAG) { if(pseudo_lookupRotamer[i_res_rot].sideAtoms[ii].atom_ptr->donorflag[0] != ' ') { pseudo_lookupRotamer[i_res_rot].sidechain_hbond_status[jj] = pseudo_lookupRotamer[i_res_rot].sideAtoms[ii].hbond_satisfied; ++jj; } ++ii; } } else pseudo_lookupRotamer[i_res_rot].sidechain_hbond_status = NULL; /* calculate internal rotamer cross-polarization energies, now that the born radii have been calculated */ E_pol_intra = 0; if(GB_FLAG==1) { VRESROT.coulombic = VRESROT.first_coulombic; while(VRESROT.coulombic->next!=NULL) { alpha = pseudo_lookupRotamer[i_res_rot].sideAtoms[VRESROT.coulombic->index_i].born_radius*pseudo_lookupRotamer[i_res_rot].sideAtoms[VRESROT.coulombic->index_j].born_radius; f = f_of_r2_alpha(VRESROT.coulombic->r2, alpha); EE = CHARGE_PRODUCT[pseudo_lookupRotamer[i_res_rot].sideAtoms[VRESROT.coulombic->index_i].atom_ptr->charge_class][pseudo_lookupRotamer[i_res_rot].sideAtoms[VRESROT.coulombic->index_j].atom_ptr->charge_class]/f; if(VRESROT.coulombic->flag == 2) { if(EE > POL_ATTRACTION_CAP) EE = POL_ATTRACTION_CAP; } E_pol_intra += BORN_CONST(f)*EE; VRESROT.coulombic = VRESROT.coulombic->next; } VRESROT.coulombic = VRESROT.first_coulombic; E_pol_intra = COULOMB_CONST*E_pol_intra; } /* local rotamer relaxation for vdw energies */ pseudo_lookupRotamer[i_res_rot].rotamerlet = NULL; calc_var_fix_energy_flag = 1; free_ligand_flag = 0; if(varPos[i].choice[i_res].resparam_ptr->ligand_flag == 1) { /* ligand rotamer that makes no interactions w/ backbone; probably a free ligand, so don't waste time and memory */ if(fabs(pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_vdw + pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_coulomb) <= EPS) if(i_res_rot != 1) { calc_var_fix_energy_flag = 0; free_ligand_flag = 1; } } /* add in intrinsic rotamer energies */ pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_1_4 += ((double)TORSION_FLAG)*VRESROT.energy.E_1_4; pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_coulomb += ((double)COULOMB_FLAG)*VRESROT.energy.E_coulomb; pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_pol += E_pol_intra; energy_rotamerlet.E_vdw = pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_vdw; if(LOCAL_MINIMIZATION_FLAG==1 && calc_var_fix_energy_flag==1 ) { pseudo_lookupRotamer[i_res_rot].rotamerlet = (ROTAMERLET *)calloc((varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->num_rotamerlets+1), sizeof(ROTAMERLET)); /* generate coords for the rotamerlets */ energy_rotamerlet.E_vdw = DBL_MAX; for(h=1;h<=varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->num_rotamerlets;++h) { build_a_SideChain(&side, &varPos[i].bkbn, varPos[i].choice[i_res].resparam_ptr, &varPos[i].seq_position, VRES->rotamerlib_ptr->rotamer[i_res_rot].rotamerlet[h]); pseudo_lookupRotamer[i_res_rot].rotamerlet[h].coord = (CARTESIAN *)calloc((VRES->numberofSideAtoms+2), sizeof(CARTESIAN)); pseudo_lookupRotamer[i_res_rot].rotamerlet[h].coord[1] = pseudo_lookupRotamer[i_res_rot].sideAtoms[1].coord; /* CB */ k=9; j=2; while(side.atom[k].atom_ptr!=NULL) { pseudo_lookupRotamer[i_res_rot].rotamerlet[h].coord[j] = side.atom[k].coord; ++k; ++j; } /* energy of the rotamerlet */ dummy_energy_rotamerlet = var_fixed_vdw_screen(fixed_atoms, pseudo_lookupRotamer[i_res_rot].rotamerlet[h].coord, VRES->atom_ptr, varPos[i].seq_position); /* identify the best rotamerlet */ if(dummy_energy_rotamerlet.E_vdw - energy_rotamerlet.E_vdw <= EPS) { best_rotamerlet = h; energy_rotamerlet = dummy_energy_rotamerlet; } /* generate coordinates for the methyl groups for this rotamerlet */ pseudo_lookupRotamer[i_res_rot].rotamerlet[h].methyl_group = NULL; if(VRES->num_methyl_groups != 0) { pseudo_lookupRotamer[i_res_rot].rotamerlet[h].methyl_group = (METHYL *)malloc(sizeof(METHYL)); pseudo_lookupRotamer[i_res_rot].rotamerlet[h].methyl_group->non_MeH = (CARTESIAN *)calloc((VRES->numberofSideAtoms+2 - 3*VRES->num_methyl_groups), sizeof(CARTESIAN)); x=1; y=1; /* identify non-Me atoms */ while(VRES->atom_ptr[x]!=NULL) { if(VRES->atom_ptr[x]->methyl_flag==0) { pseudo_lookupRotamer[i_res_rot].rotamerlet[h].methyl_group->non_MeH[y] = pseudo_lookupRotamer[i_res_rot].rotamerlet[h].coord[x]; ++y; } ++x; } pseudo_lookupRotamer[i_res_rot].rotamerlet[h].methyl_group->MeH = (METHYL_H *)calloc(VRES->num_methyl_groups+1,sizeof(METHYL_H)); for(x=1;x<=VRES->num_methyl_groups;++x) { for(y=1;y<=3;++y) { if(y==1) dummy = VRES->MeH_atom_ptr[x][1]->contactAtomNumber; else if(y==2) dummy = VRES->MeH_atom_ptr[x][1]->twobondAtomNumber; else dummy = VRES->MeH_atom_ptr[x][1]->dihedAtomNumber; if(dummy < 7) { if(dummy == 6) foundation_coords[y] = &varPos[i].bkbn.CA; else foundation_coords[y] = &varPos[i].bkbn.N; } else foundation_coords[y] = &pseudo_lookupRotamer[i_res_rot].rotamerlet[h].coord[dummy-7]; } pseudo_lookupRotamer[i_res_rot].rotamerlet[h].methyl_group->MeH[x].coord = (CARTESIAN **)calloc(4, sizeof(CARTESIAN *)); k=1; for(y=-1;y<=1;++y) { pseudo_lookupRotamer[i_res_rot].rotamerlet[h].methyl_group->MeH[x].coord[k] = (CARTESIAN *)calloc(4, sizeof(CARTESIAN)); for(j=1;j<=3;++j) { dihedangle = VRES->MeH_atom_ptr[x][j]->dihedral + METHYL_H_SPREAD*y; pseudo_lookupRotamer[i_res_rot].rotamerlet[h].methyl_group->MeH[x].coord[k][j] = chi2xyz(foundation_coords[1],foundation_coords[2],foundation_coords[3], &VRES->MeH_atom_ptr[x][j]->bondlength, &VRES->MeH_atom_ptr[x][j]->bondangle, &dihedangle); } ++k; } } } } /* tweak Methyl hydrogens for best rotamerlet if appropriate */ if(VRES->num_methyl_groups != 0) { /* energy for non-Me atoms */ energy_rotamerlet = var_fixed_vdw_screen(fixed_atoms, pseudo_lookupRotamer[i_res_rot].rotamerlet[best_rotamerlet].methyl_group->non_MeH, VRES->non_MeH_atom_ptr, varPos[i].seq_position); for(x=1;x<=VRES->num_methyl_groups;++x) { energy_Me.E_vdw = DBL_MAX; /* find best conformation for each Me */ for(y=1;y<=3;++y) { dummy_energy_Me = var_fixed_vdw_screen(fixed_atoms, pseudo_lookupRotamer[i_res_rot].rotamerlet[best_rotamerlet].methyl_group->MeH[x].coord[y], VRES->MeH_atom_ptr[x], varPos[i].seq_position); if(dummy_energy_Me.E_vdw - energy_Me.E_vdw <= EPS) { energy_Me = dummy_energy_Me; best_methyl = y; } } energy_rotamerlet.E_vdw += energy_Me.E_vdw; } } } pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_vdw = energy_rotamerlet.E_vdw; pseudo_lookupRotamer[i_res_rot].E_forcefield = energy_rotamerlet.E_vdw + VRESROT.energy.E_vdw + ((double)TORSION_FLAG)*pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_1_4 + ((double)COULOMB_FLAG)*pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_coulomb; pseudo_lookupRotamer[i_res_rot].E_gbsa = ((double)GB_FLAG)*(pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_born + pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_pol); pseudo_lookupRotamer[i_res_rot].E_asp = 0; pseudo_lookupRotamer[i_res_rot].sasa_total = 0; pseudo_lookupRotamer[i_res_rot].sasa_hydrophobic = 0; /* DUDE printf("%s\t%c\t%lf\n",varPos[i].seqpos_text_map_ptr->seqpos_text, varPos[i].choice[i_res].resparam_ptr->one_letter_code[0],pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_vdw); */ calc_var_fix_energy_flag = 0; /* anything larger is unlikely to be in the minimum */ if(pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_vdw - varPos[i].choice[i_res].resparam_ptr->max_vdw <= EPS || varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers == 1) calc_var_fix_energy_flag = 1; if(varPos[i].choice[i_res].resparam_ptr->ligand_flag == 1) { /* ligand rotamer that makes no interactions w/ backbone; probably a free ligand, so don't waste time and memory */ if(free_ligand_flag==1) { calc_var_fix_energy_flag = 0; /* if this is a free ligand, calc the energy for at least the first rotamer */ if(i_res_rot == 1) calc_var_fix_energy_flag = 1; } } // this is a position-specific rotamer if(varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot].seqpos_text!=NULL) { // this rotamer is not for this position, so don't use if(strcmp(varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot].seqpos_text, varPos[i].seqpos_text_map_ptr->seqpos_text)!=0) { calc_var_fix_energy_flag=0; pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_vdw = varPos[i].choice[i_res].resparam_ptr->max_vdw + 1e6; pseudo_lookupRotamer[i_res_rot].keep_flag = 0; } } if(varPos[i].choice[i_res].resparam_ptr->ligand_flag==1) if(varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot].native_rotamer_flag==1) { calc_var_fix_energy_flag=0; pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_vdw = varPos[i].choice[i_res].resparam_ptr->max_vdw + 1e6; pseudo_lookupRotamer[i_res_rot].keep_flag = 0; } if(calc_var_fix_energy_flag == 1) { if(SASA_FLAG!=0) /* calculate surface area for this rotamer in the empty bkbn + pseudoatoms */ { h=N_position_pseudo; while(dummyminipdb[h].seq_position == dummyminipdb[N_position_pseudo].seq_position) { dummyminipdb[h].hbond_satisfied=0; ++h; } /* attach rotamer atoms to the end of dummyminipdb, the working mini_pdbATOM array for SASA calcs */ n=num_heavy_fixedplus1; end_of_sidechain_index=num_heavy_fixedplus1; q=1; while(pseudo_lookupRotamer[i_res_rot].sideAtoms[q].seq_position!=ENDFLAG) { dummyminipdb[n]=pseudo_lookupRotamer[i_res_rot].sideAtoms[q]; dummyminipdb[n].hbond_satisfied=0; end_of_sidechain_index = n; ++n; ++q; } dummyminipdb[n].atom_ptr= NULL; dummyminipdb[n].sasa=ENDFLAG; dummyminipdb[n].seq_position=ENDFLAG; SASAcalc(dummyminipdb); /* actual SASA calculation */ exposed_hbond_groups(dummyminipdb); /* get the SASA for sidechain atoms for this rotamer */ x=1; jj=1; h=num_heavy_fixedplus1; /* num_heavy_fixedplus1 = index of first rotamer atom in dummyminipdb */ while(h<=end_of_sidechain_index) { pseudo_lookupRotamer[i_res_rot].sasa_total += dummyminipdb[h].sasa; pseudo_lookupRotamer[i_res_rot].E_asp += dummyminipdb[h].sasa*dummyminipdb[h].atom_ptr->atom_ptr->asp; if(dummyminipdb[h].atom_ptr->atom_ptr->atomtype[0] == 'C' || dummyminipdb[h].atom_ptr->atom_ptr->atomtype[0] == 'H' || dummyminipdb[h].atom_ptr->atom_ptr->atomtype[0] == 'S') pseudo_lookupRotamer[i_res_rot].sasa_hydrophobic += dummyminipdb[h].sasa; while(pseudo_lookupRotamer[i_res_rot].sideAtoms[x].atom_ptr->atomnumber != dummyminipdb[h].atom_ptr->atomnumber) ++x; pseudo_lookupRotamer[i_res_rot].sideAtoms[x].sasa = dummyminipdb[h].sasa; if(dummyminipdb[h].atom_ptr->donorflag[0] != ' ') { pseudo_lookupRotamer[i_res_rot].sidechain_hbond_status[jj] = dummyminipdb[h].hbond_satisfied || pseudo_lookupRotamer[i_res_rot].sidechain_hbond_status[jj]; ++jj; } ++h; } /* advance to the bkbn part of the pseudo_lookupRotamer[i_res_rot].sideAtoms[x] array */ while(pseudo_lookupRotamer[i_res_rot].sideAtoms[x].seq_position!=ENDFLAG) ++x; ++x; /* get the SASA for backbone atoms for this rotamer */ jj=N_position_bkbn_hbond_status; h=N_position_pseudo; /* N_position_pseudo = index of first bkbn atom for this residue in dummyminipdb */ while(dummyminipdb[h].seq_position == dummyminipdb[N_position_pseudo].seq_position) { pseudo_lookupRotamer[i_res_rot].sasa_total += dummyminipdb[h].sasa; pseudo_lookupRotamer[i_res_rot].E_asp += dummyminipdb[h].sasa*dummyminipdb[h].atom_ptr->atom_ptr->asp; if(dummyminipdb[h].atom_ptr->atom_ptr->atomtype[0] == 'C' || dummyminipdb[h].atom_ptr->atom_ptr->atomtype[0] == 'H' || dummyminipdb[h].atom_ptr->atom_ptr->atomtype[0] == 'S') pseudo_lookupRotamer[i_res_rot].sasa_hydrophobic += dummyminipdb[h].sasa; while(pseudo_lookupRotamer[i_res_rot].sideAtoms[x].atom_ptr->atomnumber != dummyminipdb[h].atom_ptr->atomnumber) ++x; pseudo_lookupRotamer[i_res_rot].sideAtoms[x].sasa = dummyminipdb[h].sasa; if(dummyminipdb[h].atom_ptr->donorflag[0] != ' ') { /* printf("%d\t%s\t%d\t%d\t%lf\t",dummyminipdb[h].seq_position, dummyminipdb[h].atom_ptr->atomname, pseudo_lookupRotamer[i_res_rot].bkbn_hbond_status[jj], dummyminipdb[h].hbond_satisfied, dummyminipdb[h].sasa); */ pseudo_lookupRotamer[i_res_rot].bkbn_hbond_status[jj] = dummyminipdb[h].hbond_satisfied || pseudo_lookupRotamer[i_res_rot].bkbn_hbond_status[jj]; /* printf("%d\n",pseudo_lookupRotamer[i_res_rot].bkbn_hbond_status[jj]); */ ++jj; } ++h; } pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_sasa = GENERAL_ASP*pseudo_lookupRotamer[i_res_rot].sasa_total; pseudo_lookupRotamer[i_res_rot].E_gbsa += pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_sasa; pseudo_lookupRotamer[i_res_rot].E_gbsa += HYDROPHOBIC_ASP*pseudo_lookupRotamer[i_res_rot].sasa_hydrophobic; } } /* DUDE printf("%d\t%c\t%lf\t%lf\n",varPos[i].seq_position, varPos[i].choice[i_res].resparam_ptr->one_letter_code[0], pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_born, HYDROPHOBIC_ASP*pseudo_lookupRotamer[i_res_rot].sasa_hydrophobic); */ pseudo_lookupRotamer[i_res_rot].keep_flag = 0; /* keep this rotamer since its vdw and solvation energy are less than the cutoff */ if(varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers == 1 || (pseudo_lookupRotamer[i_res_rot].E_gbsa - varPos[i].choice[i_res].resparam_ptr->max_solvation <= EPS && pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_vdw - varPos[i].choice[i_res].resparam_ptr->max_vdw <= EPS)) { pseudo_lookupRotamer[i_res_rot].keep_flag = 1; ++numRot; } if(i==1) /* arbitrarily attach E_sasa for pro/gly to the first variable residue */ { pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_sasa += E_sasa_pro_gly; pseudo_lookupRotamer[i_res_rot].E_asp += E_sasa_pro_gly_asp; pseudo_lookupRotamer[i_res_rot].sasa_total += sasa_pro_gly_total; pseudo_lookupRotamer[i_res_rot].sasa_hydrophobic += sasa_pro_gly_hydrophobic; pseudo_lookupRotamer[i_res_rot].E_gbsa += E_sasa_pro_gly + ((double)SASA_FLAG)*HYDROPHOBIC_ASP*sasa_pro_gly_hydrophobic; } /* identify the rotamer with the best vdw energy */ if(pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_vdw - E_forcefield_min <= EPS) { best_rotamer=i_res_rot; E_forcefield_min = pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_vdw; E_gbsa_min = pseudo_lookupRotamer[i_res_rot].E_gbsa; } } /* end i_res_rot */ /* At this point, reject rotamers that suck; copy the rest into the actual lookup table. Also calculate rotamer probabilities for this residue; use these to bias moves in MC or GA */ if(numRot==0) /* if all the rotamers for this residue suck, keep the best one */ { numRot = 1; pseudo_lookupRotamer[best_rotamer].keep_flag = 1; } else { /* calc boltz probabs based on the vdw energy for each rotamer for this residue at this position */ boltz_denom = 0; for(i_res_rot=1; i_res_rot<= varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot) { if(pseudo_lookupRotamer[i_res_rot].keep_flag == 1) { pseudo_lookupRotamer[i_res_rot].P = exp(-(pseudo_lookupRotamer[i_res_rot].energy_var_fix.E_vdw - pseudo_lookupRotamer[best_rotamer].energy_var_fix.E_vdw)); boltz_denom += pseudo_lookupRotamer[i_res_rot].P; } } for(i_res_rot=1; i_res_rot<= varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot) { if(pseudo_lookupRotamer[i_res_rot].keep_flag == 1) { /* reject rotamers that have E_vdwrotamerlib_ptr->number_of_rotamers; rotamerlib_ptr = varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr; varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr = (ROTAMERLIB *)malloc(sizeof(ROTAMERLIB)); *varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr = *rotamerlib_ptr; varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers = numRot; varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer = (ROTAMER *)calloc((numRot+2), sizeof(ROTAMER)); varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->freq_array = (double *)calloc((numRot+2), sizeof(double)); varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->freq_array[0] = 0; varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->freq_array[1] = 0; i_res_rot = 0; denom=0; for(q=1; q<= numRot0; ++q) { if(pseudo_lookupRotamer[q].keep_flag == 1) { ++i_res_rot; varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot] = *pseudo_lookupRotamer[q].rotamer; denom += varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot].freq; } } i_res_rot = 0; for(q=1; q<= numRot0; ++q) if(pseudo_lookupRotamer[q].keep_flag == 1) { ++i_res_rot; varPos[i].choice[i_res].bkbn = &varPos[i].bkbn; varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot].rank = i_res_rot; varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot].freq = varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot].freq/denom; varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->freq_array[i_res_rot] = varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->freq_array[i_res_rot-1] + varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot].freq; /* special for ligand rotamers defined by coordinates */ if(varPos[i].choice[i_res].resparam_ptr->ligand_flag==1) { if(varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->numChi==1) varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot].chi[1] = i_res_rot; } } i_res_rot = 0; for(q=1; q<= numRot0; ++q) { if(pseudo_lookupRotamer[q].keep_flag == 1) { ++i_res_rot; ERESROT.rotamer = &varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot]; ERESROT.P = NULL; ERESROT.sidechain_hbond_status = NULL; /* copy rotamer coords to milli_pdbATOM; save memory compared to mini_pdbATOM */ ERESROT.sideAtoms = (milli_pdbATOM *)calloc((varPos[i].choice[i_res].resparam_ptr->numberofSideAtoms+2), sizeof(milli_pdbATOM)); h=1; while(pseudo_lookupRotamer[q].sideAtoms[h].seq_position!=ENDFLAG) { ERESROT.sideAtoms[h].coord = pseudo_lookupRotamer[q].sideAtoms[h].coord; ERESROT.sideAtoms[h].born_radius = pseudo_lookupRotamer[q].sideAtoms[h].born_radius; /* printf("%d\t%s\t%lf\t%lf\n",varPos[i].seq_position, pseudo_lookupRotamer[q].sideAtoms[h].atom_ptr->atomname, pseudo_lookupRotamer[q].sideAtoms[h].born_radius, pseudo_lookupRotamer[q].sideAtoms[h].sasa); */ /* if(pseudo_lookupRotamer[q].sideAtoms[h].atom_ptr->donorflag[0]!=' ') printf("%s\t%lf\t%d\n",pseudo_lookupRotamer[q].sideAtoms[h].atom_ptr->atomname, pseudo_lookupRotamer[q].sideAtoms[h].sasa,pseudo_lookupRotamer[q].sideAtoms[h].hbond_satisfied); */ ++h; } ERESROT.num_var_fix_hbond=0; ERESROT.bkbn_hbond_indicies = NULL; for(x=1;x<=num_fixed_atom_hbonding_atoms;++x) { if(pseudo_lookupRotamer[q].bkbn_hbond_status[x] == 1) // this residue hbonds with or exposes bkbn if(fixed_atoms_hbond_satisfaction_array[x] == 0) // not already satisfied ++ERESROT.num_var_fix_hbond; } if(ERESROT.num_var_fix_hbond!=0) { ERESROT.bkbn_hbond_indicies = (short int *)calloc(ERESROT.num_var_fix_hbond,sizeof(short int)); z=0; for(x=1;x<=num_fixed_atom_hbonding_atoms;++x) { if(pseudo_lookupRotamer[q].bkbn_hbond_status[x] == 1) if(fixed_atoms_hbond_satisfaction_array[x] == 0) { ERESROT.bkbn_hbond_indicies[z] = x; ++z; } } } /* printf("%d\t%s\t%d\t%d\t",varPos[i].seq_position,varPos[i].choice[i_res].resparam_ptr->residuetype, count_false(pseudo_lookupRotamer[q].bkbn_hbond_status,1,num_fixed_atom_hbonding_atoms),numRot); fprintf_bool_array(stdout,pseudo_lookupRotamer[q].bkbn_hbond_status,1,num_fixed_atom_hbonding_atoms); */ free_memory(pseudo_lookupRotamer[q].bkbn_hbond_status); if(varPos[i].choice[i_res].resparam_ptr->num_hbonding_sidechain_atoms!=0) { ERESROT.sidechain_hbond_status = pseudo_lookupRotamer[q].sidechain_hbond_status; /* printf("%d\t%s\t%d\t",varPos[i].seq_position,varPos[i].choice[i_res].resparam_ptr->residuetype, count_false(ERESROT.sidechain_hbond_status,1,varPos[i].choice[i_res].resparam_ptr->num_hbonding_sidechain_atoms)); fprintf_bool_array(stdout,ERESROT.sidechain_hbond_status,1,varPos[i].choice[i_res].resparam_ptr->num_hbonding_sidechain_atoms); */ } ERESROT.rotamerlet = pseudo_lookupRotamer[q].rotamerlet; ERESROT.energy_var_fix = pseudo_lookupRotamer[q].E_forcefield + pseudo_lookupRotamer[q].E_gbsa; lookupEnergy[i].lookupRes[i_res].E_transfer += pseudo_lookupRotamer[q].P*pseudo_lookupRotamer[q].E_asp; lookupEnergy[i].lookupRes[i_res].sasa_hphob += pseudo_lookupRotamer[q].P*pseudo_lookupRotamer[q].sasa_hydrophobic; ERESROT.sasa_hphob = pseudo_lookupRotamer[q].sasa_hydrophobic; lookupEnergy[i].lookupRes[i_res].sasa_total += pseudo_lookupRotamer[q].P*pseudo_lookupRotamer[q].sasa_total; /* allocate memory for LOOKUP_ENERGY_X */ ERESROT.lookupX = (LOOKUP_ENERGY_X *)calloc((numVarPositions-i + 1), sizeof(LOOKUP_ENERGY_X)); } else { if(pseudo_lookupRotamer[q].rotamerlet!=NULL) free_ROTAMERLET(pseudo_lookupRotamer[q].rotamerlet, varPos[i].choice[i_res].resparam_ptr); if(pseudo_lookupRotamer[q].bkbn_hbond_status != NULL) free_memory(pseudo_lookupRotamer[q].bkbn_hbond_status); if(pseudo_lookupRotamer[q].sidechain_hbond_status != NULL) free_memory(pseudo_lookupRotamer[q].sidechain_hbond_status); } free_memory(pseudo_lookupRotamer[q].sideAtoms); free_memory(pseudo_lookupRotamer[q].rotamer); } save_lookupRes_to_disk(lookupEnergy[i].seq_position, varPos[i].fixed_flag, num_fixed_atom_hbonding_atoms, &lookupEnergy[i].lookupRes[i_res], varPos[i].choice[i_res].resparam_ptr, parameters); } /* end calculation for this i_res */ } /* end i_res */ n=1; q=N_position_pseudo; /* replace backbone atoms w/ pseudoatoms in pseudo_fixed_atoms */ while(holding_pseudo_fixed_atoms[n].seq_position!=ENDFLAG) { pseudo_fixed_atoms[q] = holding_pseudo_fixed_atoms[n]; ++n; ++q; } } /* end i */ /* free memory of structures used for var-fix calculations */ var_fixed_energy_calc(fixed_atoms, fixed_atoms, 2); // free the coulombic linked-list free_memory(fixed_atoms); fixed_atoms = NULL; free_memory(pseudo_sidechain); pseudo_sidechain = NULL; free_memory(dummyminipdb); dummyminipdb = NULL; free_memory(side_minipdb); side_minipdb = NULL; free_memory(holding_pseudo_fixed_atoms); holding_pseudo_fixed_atoms = NULL; free_memory(actual_CB_atom); actual_CB_atom = NULL; free_memory(foundation_coords); foundation_coords = NULL; for(i=0;iatom_ptr!=NULL) { free_memory(pseudo_fixed_atoms[i].atom_ptr->atom_ptr); pseudo_fixed_atoms[i].atom_ptr->atom_ptr = NULL; } free_memory(pseudo_fixed_atoms[i].atom_ptr); pseudo_fixed_atoms[i].atom_ptr = NULL; } free_memory(pseudo_fixed_atoms); pseudo_fixed_atoms = NULL; free_memory(pro_and_gly_pos); pro_and_gly_pos = NULL; free_memory(pseudo_lookupRotamer); pseudo_lookupRotamer = NULL; free_memory(bkbn); bkbn = NULL; chr_i.genes = chr_i.firstgene; free_memory(chr_i.genes->nextgene); chr_i.genes->nextgene = NULL; free_memory(chr_i.genes); chr_i.genes = NULL; } /* DUDE exit(0); */ /* subtract reference state energies - these are NOT stored to disk, so must be added each time */ /* also calculate final rotamer complexity */ rotamer_complexity=0; for(i=1;i<=numVarPositions; ++i) { varPos[i].dead_ended_flag=0; num_rot_position=0; for(i_res=1; i_res<=varPos[i].number_of_choices; ++i_res) { lookupEnergy[i].lookupRes[i_res].fixed_flag = varPos[i].fixed_flag; num_rot_position+=varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; varPos[i].choice[i_res].in_use_flag=1; /* default to all choices in use */ for(i_res_rot=1; i_res_rot<= varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; ++i_res_rot) { ERESROT.in_use_flag = 1; /* avoid native rotamers by penalizing their self energy */ if(protein->parameters.avoid_wt_rotamer_flag==1) if(varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->rotamer[i_res_rot].native_rotamer_flag==1) if(varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers > 1) { ERESROT.in_use_flag = 0; /* should prevent this from being selected */ ERESROT.energy_var_fix += 1e10; /* in case it doesn't */ } /* for PAIR_ENERGY_TABLE; use for when LOCAL_MINIMIZATION_FLAG == 1; for others, use JOBTYPE ENERGY_PROFILE; this is not relavant for the latter */ if(PAIR_ENERGY_TABLE_FLAG == 0) { ERESROT.E_reference = varPos[i].choice[i_res].resparam_ptr->E_reference + i_res_rot_to_specificity_energy(&varPos[i], i_res, i_res_rot); ERESROT.energy_var_fix = ERESROT.energy_var_fix - ERESROT.E_reference; } } } rotamer_complexity += log10(num_rot_position); } /* identify interacting sidechain pairs; if they don't interact, flag so memory is not reserved for them and time is not wasted calculating their energies */ if(numVarPositions > 1) { calc_var_var_flag=1; sprintf(directoryname,"interaction_lists"); while(create_directory(directoryname, parameters->lookup_energy_table_directory)==0); sprintf(directoryname,"coordinates"); while(create_directory(directoryname, parameters->lookup_energy_table_directory)==0); /* load lookupResRes from disk; contains info on who interacts with whom */ /* returns 1 if there are residue pairs for which this info is not yet available */ if(fixed_positions_flag == 0) calc_var_var_flag = load_lookupResRes_from_disk(lookupEnergy, varPos, parameters); if(calc_var_var_flag==1) { i_atoms = (milli_pdbATOM **)calloc(MAX_ROTAMERS, sizeof(milli_pdbATOM *)); j_atoms = (milli_pdbATOM **)calloc(MAX_ROTAMERS, sizeof(milli_pdbATOM *)); for(i=1;i<=numVarPositions; ++i) { xmax=varPos[i].bkbn.CB.x+FORCEFIELD_DISTANCE_CUTOFF; xmin=varPos[i].bkbn.CB.x-FORCEFIELD_DISTANCE_CUTOFF; ymax=varPos[i].bkbn.CB.y+FORCEFIELD_DISTANCE_CUTOFF; ymin=varPos[i].bkbn.CB.y-FORCEFIELD_DISTANCE_CUTOFF; zmax=varPos[i].bkbn.CB.z+FORCEFIELD_DISTANCE_CUTOFF; zmin=varPos[i].bkbn.CB.z-FORCEFIELD_DISTANCE_CUTOFF; for(i_res=1; i_res<=varPos[i].number_of_choices; ++i_res) { if(lookupEnergy[i].lookupRes[i_res].lookupResRes == NULL) { lookupEnergy[i].lookupRes[i_res].lookupResRes = (LOOKUP_ENERGY_RESIDUE_RESIDUE *)calloc((numVarPositions-i + 1), sizeof(LOOKUP_ENERGY_RESIDUE_RESIDUE)); for(j=(i+1);j<=numVarPositions;++j) { lookupEnergy[i].lookupRes[i_res].lookupResRes[j-i].resres_interact_flag = NULL; /* do CB's interact? if so, then all the sidechains for these positions interact */ if(varPos[i].fixed_flag==0 || varPos[j].fixed_flag==0) /* at least one of the positions must be floating */ if(varPos[j].bkbn.CB.x<=xmax) if(varPos[j].bkbn.CB.x>=xmin) if(varPos[j].bkbn.CB.y<=ymax) if(varPos[j].bkbn.CB.y>=ymin) if(varPos[j].bkbn.CB.z<=zmax) if(varPos[j].bkbn.CB.z>=zmin) if(Distance_sqrd(varPos[j].bkbn.CB,varPos[i].bkbn.CB) <= FORCEFIELD_DISTANCE_CUTOFF_SQRD) { lookupEnergy[i].lookupRes[i_res].lookupResRes[j-i].resres_interact_flag = (int *)calloc(varPos[j].number_of_choices, sizeof(int)); for(h=0;hrotamerlib_ptr->number_of_rotamers; ++i_res_rot) { if(ERESROT.sideAtoms == NULL) { if(varPos[i].fixed_flag == 0) sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure", parameters->lookup_energy_table_directory, lookupEnergy[i].seq_position, varPos[i].choice[i_res].resparam_ptr->residuetype, lookupEnergy[i].seq_position); else sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure.fixed", parameters->lookup_energy_table_directory, lookupEnergy[i].seq_position, varPos[i].choice[i_res].resparam_ptr->residuetype, lookupEnergy[i].seq_position); load_coordinates_lookupRes(&lookupEnergy[i].lookupRes[i_res], varPos[i].choice[i_res].resparam_ptr, dummystring); } i_atoms[i_res_rot] = ERESROT.sideAtoms; } i_atoms[varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers + 1] = NULL; for(j=(i+1);j<=numVarPositions;++j) { if(varPos[i].fixed_flag==0 || varPos[j].fixed_flag==0) /* at least one of the positions must be floating */ if(lookupEnergy[i].lookupRes[i_res].lookupResRes[j-i].resres_interact_flag == NULL) /* if not NULL, CB found to interact above */ { lookupEnergy[i].lookupRes[i_res].lookupResRes[j-i].resres_interact_flag = (int *)calloc(varPos[j].number_of_choices, sizeof(int)); for(h=0;hrotamerlib_ptr->number_of_rotamers; ++j_res_rot) { if(ERESROTj.sideAtoms == NULL) { if(varPos[j].fixed_flag == 0) sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure", parameters->lookup_energy_table_directory, lookupEnergy[j].seq_position, varPos[j].choice[j_res].resparam_ptr->residuetype, lookupEnergy[j].seq_position); else sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure.fixed", parameters->lookup_energy_table_directory, lookupEnergy[j].seq_position, varPos[j].choice[j_res].resparam_ptr->residuetype, lookupEnergy[j].seq_position); load_coordinates_lookupRes(&lookupEnergy[j].lookupRes[j_res], varPos[j].choice[j_res].resparam_ptr, dummystring); } j_atoms[j_res_rot] = ERESROTj.sideAtoms; } j_atoms[varPos[j].choice[j_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers + 1] = NULL; /* interaction */ if(do_these_interact(i_atoms, j_atoms, (varPos[i].choice[i_res].resparam_ptr->numberofSideAtoms+1), (varPos[j].choice[j_res].resparam_ptr->numberofSideAtoms+1)) == 1) { lookupEnergy[i].lookupRes[i_res].lookupResRes[j-i].resres_interact_flag[j_res-1] = 1; } } } } if(fixed_positions_flag == 0) /* save this data to disk */ save_lookupResRes_to_disk(i, i_res, lookupEnergy, varPos, parameters); } } } free_memory(i_atoms); free_memory(j_atoms); i_atoms = NULL; j_atoms = NULL; } /* if(parameters->disk_lookup_table_flag == 1) // generate the whole pair table at once { load_lookupResX_from_disk(lookupEnergy, varPos, parameters); } else // generate the pair table as needed; for rotamer_calc_foremen; set up the table but don't fill it in */ load_lookupResX_from_disk(lookupEnergy, varPos, parameters); // if parameters->disk_lookup_table_flag=1, pair table is generated completely; otherwise, need to make space if(parameters->disk_lookup_table_flag != 1) for(i=1;i<=numVarPositions; ++i) /* position i loop */ { for(i_res=1; i_res<=varPos[i].number_of_choices; ++i_res) /* residuetype i_res loop */ { num_rot_position = num_rot_position+lookupEnergy[i].lookupRes[i_res].sizeof_lookupRot; for(i_res_rot=1; i_res_rot<= lookupEnergy[i].lookupRes[i_res].sizeof_lookupRot; ++i_res_rot) { /* position j loop */ for(j=(i+1);j<=numVarPositions;++j) { go_flag = 0; if(!(varPos[i].fixed_flag == 1 || varPos[j].fixed_flag == 1)) go_flag = 1; if(go_flag==1) { if(lookupEnergy[i].lookupRes[i_res].lookupResRes[j-i].resres_interact_flag != NULL) /* these positions interact */ { /* allocate memory for LOOKUP_ENERGY_RESIDUE_X */ ERESROT.lookupX[j-i].lookupResX = (LOOKUP_ENERGY_RESIDUE_X *)calloc((varPos[j].number_of_choices),sizeof(LOOKUP_ENERGY_RESIDUE_X)); /* residue j_res loop */ for(j_res=0;j_res1) { i_res=1; i_res_rot=1; if(ERESROT.sideAtoms == NULL) { sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure", parameters->lookup_energy_table_directory, lookupEnergy[i].seq_position, varPos[i].choice[i_res].resparam_ptr->residuetype, lookupEnergy[i].seq_position); load_coordinates_lookupRes(&lookupEnergy[i].lookupRes[i_res], varPos[i].choice[i_res].resparam_ptr, dummystring); } varPos[i].fixed_flag=2; // not really "fixed", but acts as if it is fixed_positions_flag=1; } else ++dummy; } } } } new_num_fixed_atom_hbonding_atoms = num_fixed_atom_hbonding_atoms; FIXED_ATOM_HBOND_STATUS_ARRAYLENGTH = new_num_fixed_atom_hbonding_atoms; BKBN_HBOND_STATUS_ARRAYLENGTH = num_fixed_atom_hbonding_atoms; FIXED_HBOND_SATISF_ARRAY = (bool *)calloc(new_num_fixed_atom_hbonding_atoms+2,sizeof(bool)); falsify_bool_array(FIXED_HBOND_SATISF_ARRAY,1,new_num_fixed_atom_hbonding_atoms); for(x=1;x<=num_fixed_atom_hbonding_atoms;++x) FIXED_HBOND_SATISF_ARRAY[x] = fixed_atoms_hbond_satisfaction_array[x]; free_memory(fixed_atoms_hbond_satisfaction_array); // compress FIXED_HBOND_SATISF_ARRAY; get rid of satisfied entries for(i=1;i<=numVarPositions; ++i) for(i_res=1;i_res<=varPos[i].number_of_choices; ++i_res) for(i_res_rot=1; i_res_rot<= lookupEnergy[i].lookupRes[i_res].sizeof_lookupRot; ++i_res_rot) if(ERESROT.num_var_fix_hbond!=0) { x=0; k=ERESROT.num_var_fix_hbond; while(x < k) { if(FIXED_HBOND_SATISF_ARRAY[ERESROT.bkbn_hbond_indicies[x]] == 1) { ERESROT.bkbn_hbond_indicies[x] = 0; --ERESROT.num_var_fix_hbond; } ++x; } if(k!=ERESROT.num_var_fix_hbond) // resize the ERESROT.bkbn_hbond_indicies array { if(ERESROT.num_var_fix_hbond!=0) { holding_shortint_array = ERESROT.bkbn_hbond_indicies; ERESROT.bkbn_hbond_indicies = (short int *)calloc(ERESROT.num_var_fix_hbond,sizeof(short int)); x=0; for(j=0;jparameters.algorithm, "MC_PAIR_ENERGY_TABLE")!=0 && strstr(protein->parameters.algorithm, "SCMF")==0) { // place backbone due to fixed sidechains hbond status in FIXED_HBOND_SATISF_ARRAY for(i=1;i<=numVarPositions; ++i) if(varPos[i].fixed_flag != 0) { i_res=1; i_res_rot=1; if(ERESROT.num_var_fix_hbond!=0) { for(x=0;x < ERESROT.num_var_fix_hbond;++x) FIXED_HBOND_SATISF_ARRAY[ERESROT.bkbn_hbond_indicies[x]]=1; free_memory(ERESROT.bkbn_hbond_indicies); ERESROT.num_var_fix_hbond = 0; } } // adjust ERESROT.bkbn_hbond_indicies to match up to new indicies of affected atoms // compress FIXED_HBOND_SATISF_ARRAY; get rid of satisfied entries for(i=1;i<=numVarPositions; ++i) if(varPos[i].fixed_flag == 0) for(i_res=1;i_res<=varPos[i].number_of_choices; ++i_res) for(i_res_rot=1; i_res_rot<= lookupEnergy[i].lookupRes[i_res].sizeof_lookupRot; ++i_res_rot) if(ERESROT.num_var_fix_hbond!=0) { x=0; k=ERESROT.num_var_fix_hbond; while(x < k) { if(FIXED_HBOND_SATISF_ARRAY[ERESROT.bkbn_hbond_indicies[x]] == 1) { ERESROT.bkbn_hbond_indicies[x] = 0; --ERESROT.num_var_fix_hbond; } ++x; } if(k!=ERESROT.num_var_fix_hbond) // resize the ERESROT.bkbn_hbond_indicies array { if(ERESROT.num_var_fix_hbond!=0) { holding_shortint_array = ERESROT.bkbn_hbond_indicies; ERESROT.bkbn_hbond_indicies = (short int *)calloc(ERESROT.num_var_fix_hbond,sizeof(short int)); x=0; for(j=0;jnum_hbonding_sidechain_atoms; holding_fixed_hbond_satisf_array = FIXED_HBOND_SATISF_ARRAY; FIXED_HBOND_SATISF_ARRAY = (bool *)calloc(new_num_fixed_atom_hbonding_atoms+2,sizeof(bool)); holding_atomresparam = FIXED_ATOM_HBOND_ATOMRESPARAM; FIXED_ATOM_HBOND_ATOMRESPARAM = (ATOMRESPARAM **)calloc(new_num_fixed_atom_hbonding_atoms+2,sizeof(ATOMRESPARAM *)); j=0; for(i=1;i<=num_fixed_atom_hbonding_atoms;++i) if(holding_fixed_hbond_satisf_array[i] == 0) { ++j; FIXED_HBOND_SATISF_ARRAY[j] = holding_fixed_hbond_satisf_array[i]; FIXED_ATOM_HBOND_ATOMRESPARAM[j] = holding_atomresparam[i]; } free_memory(holding_fixed_hbond_satisf_array); free_memory(holding_atomresparam); num_fixed_atom_hbonding_atoms = j; BKBN_HBOND_STATUS_ARRAYLENGTH = num_fixed_atom_hbonding_atoms; new_num_fixed_atom_hbonding_atoms = num_fixed_atom_hbonding_atoms; for(i=1;i<=numVarPositions; ++i) if(varPos[i].fixed_flag != 0) new_num_fixed_atom_hbonding_atoms += varPos[i].choice[1].resparam_ptr->num_hbonding_sidechain_atoms; FIXED_ATOM_HBOND_STATUS_ARRAYLENGTH = new_num_fixed_atom_hbonding_atoms; // append fixed due to backbone hbond status to FIXED_HBOND_SATISF_ARRAY ii=num_fixed_atom_hbonding_atoms+1; for(i=1;i<=numVarPositions; ++i) if(varPos[i].fixed_flag != 0) { if(varPos[i].choice[1].resparam_ptr->num_hbonding_sidechain_atoms!=0) { i_res=1; i_res_rot=1; copy_bool_array(FIXED_HBOND_SATISF_ARRAY, ii, ERESROT.sidechain_hbond_status, 1, varPos[i].choice[1].resparam_ptr->num_hbonding_sidechain_atoms); j=8; x=ii; while(varPos[i].choice[1].resparam_ptr->atom[j].atomnumber!=ENDFLAG) { if(varPos[i].choice[1].resparam_ptr->atom[j].other_info>0) if(varPos[i].choice[1].resparam_ptr->atom[j].donorflag[0] !=' ') { FIXED_ATOM_HBOND_ATOMRESPARAM[x] = &varPos[i].choice[1].resparam_ptr->atom[j]; ++x; } ++j; } ii+=varPos[i].choice[1].resparam_ptr->num_hbonding_sidechain_atoms; } } for(i=1;i<=numVarPositions; ++i) { if(varPos[i].fixed_flag == 0) // i floats { for(i_res=1; i_res<=varPos[i].number_of_choices; ++i_res) { for(i_res_rot=1; i_res_rot<= lookupEnergy[i].lookupRes[i_res].sizeof_lookupRot; ++i_res_rot) { ii=num_fixed_atom_hbonding_atoms+1; for(j=1;j<=numVarPositions;++j) { if(j!=i) if(varPos[j].fixed_flag != 0) // j is fixed { j_res=1; j_res_rot=1; if(j>i) { if(ERESROT.lookupX[j-i].lookupResX != NON_INTERACT_LOOKUP_RES_X) { if(ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX != NON_INTERACT_LOOKUP_ROT_X) { if(ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX == NULL) { j_gene.choice_ptr = &varPos[j].choice[j_res]; j_gene.seq_position = varPos[j].seq_position; j_gene.j_choice_index = j_res-1; j_gene.lookupRot = &j_gene.choice_ptr->lookup_res_ptr->lookupRot[j_res_rot]; j_gene.lookupRot_index = j_res_rot; j_gene.varpos_ptr = &varPos[j]; i_gene.choice_ptr = &varPos[i].choice[i_res]; i_gene.seq_position = varPos[i].seq_position; i_gene.j_choice_index = i_res-1; i_gene.varpos_ptr = &varPos[i]; i_gene.lookupRot = &ERESROT; i_gene.lookupRot_index = i_res_rot; if(ERESROT.sideAtoms == NULL) { sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure", parameters->lookup_energy_table_directory, lookupEnergy[i].seq_position, varPos[i].choice[i_res].resparam_ptr->residuetype, lookupEnergy[i].seq_position); load_coordinates_lookupRes(&lookupEnergy[i].lookupRes[i_res], varPos[i].choice[i_res].resparam_ptr, dummystring); } if(ERESROTj.sideAtoms == NULL) { if(varPos[j].fixed_flag == 1) sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure.fixed", parameters->lookup_energy_table_directory, lookupEnergy[j].seq_position, varPos[j].choice[j_res].resparam_ptr->residuetype, lookupEnergy[j].seq_position); else sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure", parameters->lookup_energy_table_directory, lookupEnergy[j].seq_position, varPos[j].choice[j_res].resparam_ptr->residuetype, lookupEnergy[j].seq_position); load_coordinates_lookupRes(&lookupEnergy[j].lookupRes[j_res], varPos[j].choice[j_res].resparam_ptr, dummystring); } for(q=1;q<=lookupEnergy[i].lookupRes[i_res].sizeof_lookupRot;++q) { i_gene.lookupRot = &i_gene.choice_ptr->lookup_res_ptr->lookupRot[q]; i_gene.lookupRot_index = q; i_gene.lookupRot->lookupX[j-i].lookupResX[j_res-1].lookupRotX = (LOOKUP_ENERGY_ROTAMER_X *)calloc(lookupEnergy[j].lookupRes[j_res].sizeof_lookupRot,sizeof(LOOKUP_ENERGY_ROTAMER_X)); E_temp = rotamer_pair_energy(&i_gene, &j_gene, &i_gene.lookupRot->lookupX[j-i].lookupResX[j_res-1].lookupRotX_hbond ); if(E_temp == 0) i_gene.lookupRot->lookupX[j-i].lookupResX[j_res-1].lookupRotX[0].energy_var_var = NON_INTERACT_PTR; else i_gene.lookupRot->lookupX[j-i].lookupResX[j_res-1].lookupRotX[0].energy_var_var = E_temp; } trim_LOOKUP_ROTAMER_X_HBOND(&i_gene.lookupRot->lookupX[j-i].lookupResX[j_res-1].lookupRotX_hbond); } if(ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX[0].energy_var_var != NON_INTERACT_PTR) { ERESROT.energy_var_fix += ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX[j_res_rot-1].energy_var_var; if(ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX_hbond != NULL) { h=find_j_res_rot_minus1_in_hbond_rot_pair(j_res_rot-1, ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX_hbond); if(h!=ENDFLAG) { union_bool_arrays(ERESROT.sidechain_hbond_status,1, ERESROT.sidechain_hbond_status,1,varPos[i].choice[i_res].resparam_ptr->num_hbonding_sidechain_atoms, ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX_hbond->hbond_rot_pair[h].i_hbond_status, 0, varPos[i].choice[i_res].resparam_ptr->numhbondsideatoms_minus1); // add indicies of fixed sidechain atoms j that hbond with i to ERESROT.bkbn_hbond_indicies z=ii; for(x=0;xnum_hbonding_sidechain_atoms;++x) { if(ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX_hbond->hbond_rot_pair[h].j_hbond_status[x]==1) { ++ERESROT.num_var_fix_hbond; if(ERESROT.bkbn_hbond_indicies!=NULL) // ERESROT.bkbn_hbond_indicies already exists { holding_shortint_array = ERESROT.bkbn_hbond_indicies; ERESROT.bkbn_hbond_indicies = (short int *)calloc(ERESROT.num_var_fix_hbond,sizeof(short int)); y=0; while(ylookup_res_ptr->lookupRot[j_res_rot]; j_gene.lookupRot_index = j_res_rot; j_gene.varpos_ptr = &varPos[j]; i_gene.choice_ptr = &varPos[i].choice[i_res]; i_gene.seq_position = varPos[i].seq_position; i_gene.j_choice_index = i_res-1; i_gene.varpos_ptr = &varPos[i]; i_gene.lookupRot = &ERESROT; i_gene.lookupRot_index = i_res_rot; if(ERESROT.sideAtoms == NULL) { sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure", parameters->lookup_energy_table_directory, lookupEnergy[i].seq_position, varPos[i].choice[i_res].resparam_ptr->residuetype, lookupEnergy[i].seq_position); load_coordinates_lookupRes(&lookupEnergy[i].lookupRes[i_res], varPos[i].choice[i_res].resparam_ptr, dummystring); } if(ERESROTj.sideAtoms == NULL) { if(varPos[j].fixed_flag==1) sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure.fixed", parameters->lookup_energy_table_directory, lookupEnergy[j].seq_position, varPos[j].choice[j_res].resparam_ptr->residuetype, lookupEnergy[j].seq_position); else sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure", parameters->lookup_energy_table_directory, lookupEnergy[j].seq_position, varPos[j].choice[j_res].resparam_ptr->residuetype, lookupEnergy[j].seq_position); load_coordinates_lookupRes(&lookupEnergy[j].lookupRes[j_res], varPos[j].choice[j_res].resparam_ptr, dummystring); } for(q=1;q<=lookupEnergy[i].lookupRes[i_res].sizeof_lookupRot;++q) { i_gene.lookupRot = &i_gene.choice_ptr->lookup_res_ptr->lookupRot[q]; i_gene.lookupRot_index = q; E_temp = rotamer_pair_energy(&j_gene, &i_gene, &ERESROTj.lookupX[i-j].lookupResX[i_res-1].lookupRotX_hbond); if(E_temp == 0) ERESROTj.lookupX[i-j].lookupResX[i_res-1].lookupRotX[q-1].energy_var_var = NON_INTERACT_PTR; else ERESROTj.lookupX[i-j].lookupResX[i_res-1].lookupRotX[q-1].energy_var_var = E_temp; } trim_LOOKUP_ROTAMER_X_HBOND(&ERESROTj.lookupX[i-j].lookupResX[i_res-1].lookupRotX_hbond); } if(ERESROTj.lookupX[i-j].lookupResX[i_res-1].lookupRotX[i_res_rot-1].energy_var_var != NON_INTERACT_PTR) { ERESROT.energy_var_fix += ERESROTj.lookupX[i-j].lookupResX[i_res-1].lookupRotX[i_res_rot-1].energy_var_var; if(ERESROTj.lookupX[i-j].lookupResX[i_res-1].lookupRotX_hbond != NULL) { h=find_j_res_rot_minus1_in_hbond_rot_pair(i_res_rot-1, ERESROTj.lookupX[i-j].lookupResX[i_res-1].lookupRotX_hbond); if(h!=ENDFLAG) { union_bool_arrays(ERESROT.sidechain_hbond_status,1, ERESROT.sidechain_hbond_status,1,varPos[i].choice[i_res].resparam_ptr->num_hbonding_sidechain_atoms, ERESROTj.lookupX[i-j].lookupResX[i_res-1].lookupRotX_hbond->hbond_rot_pair[h].j_hbond_status, 0, varPos[i].choice[i_res].resparam_ptr->numhbondsideatoms_minus1); // add indicies of fixed sidechain atoms j that hbond with i to ERESROT.bkbn_hbond_indicies z=ii; for(x=0;xnum_hbonding_sidechain_atoms;++x) { if(ERESROTj.lookupX[i-j].lookupResX[i_res-1].lookupRotX_hbond->hbond_rot_pair[h].i_hbond_status[x]==1) { ++ERESROT.num_var_fix_hbond; if(ERESROT.bkbn_hbond_indicies!=NULL) // ERESROT.bkbn_hbond_indicies already exists { holding_shortint_array = ERESROT.bkbn_hbond_indicies; ERESROT.bkbn_hbond_indicies = (short int *)calloc(ERESROT.num_var_fix_hbond,sizeof(short int)); y=0; while(ynum_hbonding_sidechain_atoms!=0) ii += varPos[j].choice[1].resparam_ptr->num_hbonding_sidechain_atoms; } // end j } // end i_res_rot } // end i_res } else // i is fixed { i_res=1; i_res_rot=1; protein->parameters.fixedatoms_energy += ERESROT.energy_var_fix; // add i fixed energy to fixed ERESROT.energy_var_fix=FIXED_POSITION_PTR; // since the energy is now included in fixedatoms_energy if(ERESROT.sideAtoms == NULL) { if(varPos[i].fixed_flag==1) sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure.fixed", parameters->lookup_energy_table_directory, lookupEnergy[i].seq_position, varPos[i].choice[i_res].resparam_ptr->residuetype, lookupEnergy[i].seq_position); else sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure", parameters->lookup_energy_table_directory, lookupEnergy[i].seq_position, varPos[i].choice[i_res].resparam_ptr->residuetype, lookupEnergy[i].seq_position); load_coordinates_lookupRes(&lookupEnergy[i].lookupRes[i_res], varPos[i].choice[i_res].resparam_ptr, dummystring); } i_gene.nextgene = NULL; i_gene.choice_ptr = &varPos[i].choice[i_res]; i_gene.seq_position = varPos[i].seq_position; i_gene.j_choice_index = i_res-1; i_gene.lookupRot = &ERESROT; i_gene.lookupRot_index = i_res_rot; for(j=(i+1);j<=numVarPositions;++j) { if(varPos[j].fixed_flag != 0) // j and i are both fixed { j_res=1; j_res_rot=1; if(ERESROTj.sideAtoms == NULL) { if(varPos[j].fixed_flag == 1) sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure.fixed", parameters->lookup_energy_table_directory, lookupEnergy[j].seq_position, varPos[j].choice[j_res].resparam_ptr->residuetype, lookupEnergy[j].seq_position); else sprintf(dummystring, "%s/coordinates/%d/%s.%d.structure", parameters->lookup_energy_table_directory, lookupEnergy[j].seq_position, varPos[j].choice[j_res].resparam_ptr->residuetype, lookupEnergy[j].seq_position); load_coordinates_lookupRes(&lookupEnergy[j].lookupRes[j_res], varPos[j].choice[j_res].resparam_ptr, dummystring); } j_gene.nextgene = NULL; j_gene.choice_ptr = &varPos[j].choice[j_res]; j_gene.seq_position = varPos[j].seq_position; j_gene.j_choice_index = j_res-1; j_gene.lookupRot = &ERESROTj; j_gene.lookupRot_index = j_res_rot; lookupRotX_hbond=NULL; protein->parameters.fixedatoms_energy += rotamer_pair_energy(&i_gene, &j_gene, &lookupRotX_hbond); if(lookupRotX_hbond!=NULL) { union_bool_arrays(ERESROT.sidechain_hbond_status,1, ERESROT.sidechain_hbond_status,1,i_gene.choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms, lookupRotX_hbond->hbond_rot_pair[0].i_hbond_status, 0, i_gene.choice_ptr->resparam_ptr->numhbondsideatoms_minus1); union_bool_arrays(ERESROTj.sidechain_hbond_status,1, ERESROTj.sidechain_hbond_status,1,j_gene.choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms, lookupRotX_hbond->hbond_rot_pair[0].j_hbond_status, 0, j_gene.choice_ptr->resparam_ptr->numhbondsideatoms_minus1); free_LOOKUP_ROTAMER_X_HBOND(&lookupRotX_hbond); } } } free_memory(ERESROT.sideAtoms); ERESROT.sideAtoms = NULL; } } // update fixed due to backbone hbond status to FIXED_HBOND_SATISF_ARRAY ii=num_fixed_atom_hbonding_atoms+1; for(i=1;i<=numVarPositions; ++i) if(varPos[i].fixed_flag != 0) { if(varPos[i].choice[1].resparam_ptr->num_hbonding_sidechain_atoms!=0) { i_res=1; i_res_rot=1; copy_bool_array(FIXED_HBOND_SATISF_ARRAY, ii, ERESROT.sidechain_hbond_status, 1, varPos[i].choice[1].resparam_ptr->num_hbonding_sidechain_atoms); ii+=varPos[i].choice[1].resparam_ptr->num_hbonding_sidechain_atoms; free_memory(ERESROT.sidechain_hbond_status); } } /* free up memory that is no longer required */ for(i=1;i<=numVarPositions; ++i) { if(varPos[i].fixed_flag == 0) /* i floats */ { for(i_res=1; i_res<=varPos[i].number_of_choices; ++i_res) { for(i_res_rot=1; i_res_rot<= lookupEnergy[i].lookupRes[i_res].sizeof_lookupRot; ++i_res_rot) { for(j=(i+1);j<=numVarPositions;++j) { if(varPos[j].fixed_flag != 0) /* j is fixed, so its energies with i are included in i's var_fix */ { j_res=1; if(ERESROT.lookupX[j-i].lookupResX != NON_INTERACT_LOOKUP_RES_X) { if(ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX != NON_INTERACT_LOOKUP_ROT_X) if(ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX != NULL) { free_memory(ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX); ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX=NON_INTERACT_LOOKUP_ROT_X; if(ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX_hbond != NULL) free_LOOKUP_ROTAMER_X_HBOND(&ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX_hbond); } free_memory(ERESROT.lookupX[j-i].lookupResX); ERESROT.lookupX[j-i].lookupResX=NON_INTERACT_LOOKUP_RES_X; } } } } } } else /* i is fixed, so clear its entire var-var tree since the energies are already included elsewhere */ { i_res=1; i_res_rot=1; for(j=(i+1);j<=numVarPositions;++j) { if(ERESROT.lookupX[j-i].lookupResX != NON_INTERACT_LOOKUP_RES_X) { for(j_res=1;j_res<=varPos[j].number_of_choices;++j_res) { if(ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX != NULL) if(ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX != NON_INTERACT_LOOKUP_ROT_X) { free_memory(ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX); ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX=NON_INTERACT_LOOKUP_ROT_X; if(ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX_hbond != NULL) free_LOOKUP_ROTAMER_X_HBOND(&ERESROT.lookupX[j-i].lookupResX[j_res-1].lookupRotX_hbond); } } free_memory(ERESROT.lookupX[j-i].lookupResX); ERESROT.lookupX[j-i].lookupResX=NON_INTERACT_LOOKUP_RES_X; } } } } } } protein->parameters.numMovingPositions = 0; for(i=1;i<=numVarPositions; ++i) { varPos[i].number_of_resimers = 0; for(i_res=1; i_res<=varPos[i].number_of_choices; ++i_res) varPos[i].number_of_resimers += varPos[i].choice[i_res].resparam_ptr->rotamerlib_ptr->number_of_rotamers; if(varPos[i].fixed_flag == 0) ++protein->parameters.numMovingPositions; } for(i=1;i<=numVarPositions; ++i) if(varPos[i].fixed_flag == 0) /* i floats */ for(i_res=1; i_res<=varPos[i].number_of_choices; ++i_res) for(i_res_rot=1; i_res_rot<= lookupEnergy[i].lookupRes[i_res].sizeof_lookupRot; ++i_res_rot) { if(ERESROT.sideAtoms!=NULL) free_memory(ERESROT.sideAtoms); if(ERESROT.rotamerlet!=NULL) free_ROTAMERLET(ERESROT.rotamerlet, varPos[i].choice[i_res].resparam_ptr); } if(LOGFILE_FLAG!=0) { sprintf(dummystring,"%s/%d.lookup_log",CURRENT_WORKING_DIRECTORY, GET_PID); rm_file(dummystring); sprintf(dummystring,"%s.log",protein->parameters.output_prefix); now = time(NULL); logfile_ptr=fopen(dummystring,"a"); fprintf(logfile_ptr,"FILTERED_ROTAMER_COMPLEXITY\t10^%lf\t%s",rotamer_complexity,ctime(&now)); fclose(logfile_ptr); } protein->parameters.log10_rotamer_combinations=rotamer_complexity; if(MAX_MUTATIONS > 0) { chr = (CHROMOSOME *)malloc(sizeof(CHROMOSOME)); chr->genes = (MENDEL *)malloc(sizeof(MENDEL)); chr->firstgene = chr->genes; chr->bkbngenes = NULL; chr->first_bkbngene = NULL; j=1; while(varPos[j].seq_position!=ENDFLAG) { chr->genes->seq_position = varPos[j].seq_position; /* assign first residuetype */ chr->genes->varpos_ptr = &varPos[j]; chr->genes->j_choice_index = 1; CHOICE_to_GENE(chr->genes, chr->genes->varpos_ptr->choice[1], chr->genes->j_choice_index); chr->genes->nextgene = (MENDEL *)malloc(sizeof(MENDEL)); chr->genes = chr->genes->nextgene; ++j; } chr->genes->seq_position = ENDFLAG; chr->genes->nextgene=NULL; /* convert to template sequence */ sequence_to_CHROMOSOME(chr, protein->invar_pos, protein->template_sequence); generate_custom_sequence_CHROMOSOME(chr,1); // intializes stuff chr->genes = chr->firstgene; free_CHROMOSOME(chr); free_memory(chr); } free_memory(dummystring); free_memory(directoryname); free_memory(side.atom); }