/* EGAD: pairwise_energy_calc.cpp Navin Pokala and Tracy Handel Dept. of Molecular and Cell Biology University of California, Berkeley Copyright (C) 2003 Regents of the University of California GNU Public License Aug 12 2003 Absolutely no warranties are made or are implied with the use of this program or its parts. This file contains functions for calculating energies required for pair energy lookup table entries. */ #include "pairwise_energy_calc.h" /* This function returns energy for the fixed atoms. Called by input_VARIABLE_POSITIONS */ ENERGY fixedatoms_energy(mini_pdbATOM *input_fixed_atoms, VARIABLE_POSITION *varpos, RESPARAM *resparam) { int i, j, k; mini_pdbATOM *fixed_atoms; ENERGY energy; extern int GB_FLAG, COULOMB_FLAG, TORSION_FLAG; fixed_atoms = (mini_pdbATOM *)calloc(MAX_ATOMS, sizeof(mini_pdbATOM)); i=1; j=1; k=1; while(varpos[k].seq_position!=ENDFLAG) { while(input_fixed_atoms[i].seq_position!=varpos[k].seq_position && input_fixed_atoms[i].seq_position!=ENDFLAG) { fixed_atoms[j] = input_fixed_atoms[i]; ++j; ++i; } while(varpos[k].seq_position==input_fixed_atoms[i].seq_position && input_fixed_atoms[i].seq_position!=ENDFLAG) { if(input_fixed_atoms[i].atom_ptr->other_info!=-5 || strcmp(input_fixed_atoms[i].atom_ptr->atomname, "2HA")==0) { fixed_atoms[j] = input_fixed_atoms[i]; /* fixed_atoms in this function doesn't include CB */ ++j; } ++i; } ++k; } while(input_fixed_atoms[i].seq_position!=ENDFLAG) { fixed_atoms[j] = input_fixed_atoms[i]; ++j; ++i; } fixed_atoms[j].sasa=ENDFLAG; fixed_atoms[j].seq_position=ENDFLAG; /* calculate energy values and mark bkbn atoms making intra-bkbn hbonds */ energy = energy_calc(fixed_atoms, 2, 0); i=1; j=1; k=1; while(varpos[k].seq_position!=ENDFLAG) { while(input_fixed_atoms[i].seq_position!=varpos[k].seq_position && input_fixed_atoms[i].seq_position!=ENDFLAG) { input_fixed_atoms[i].hbond_satisfied = fixed_atoms[j].hbond_satisfied; ++j; ++i; } while(varpos[k].seq_position==input_fixed_atoms[i].seq_position && input_fixed_atoms[i].seq_position!=ENDFLAG) { if(input_fixed_atoms[i].atom_ptr->other_info!=-5 || strcmp(input_fixed_atoms[i].atom_ptr->atomname, "2HA")==0) { input_fixed_atoms[i].hbond_satisfied = fixed_atoms[j].hbond_satisfied; ++j; } ++i; } ++k; } while(input_fixed_atoms[i].seq_position!=ENDFLAG) { input_fixed_atoms[i].hbond_satisfied = fixed_atoms[j].hbond_satisfied; ++j; ++i; } energy.E_reference = 0; energy.num_unsatisfied_hbond = 0; free_memory(fixed_atoms); energy.E_sasa = 0; energy.E_born = ((double)GB_FLAG)*energy.E_born; energy.E_pol = ((double)GB_FLAG)*energy.E_pol; energy.E_coulomb = ((double)COULOMB_FLAG)*energy.E_coulomb; energy.E_1_4 = ((double)TORSION_FLAG)*energy.E_1_4; energy.E_hbond = ((double)HBOND_FUNCTION_FLAG)*energy.E_hbond; return(energy); } /* returns var_var_ENERGY between sidechain rotamers i and j; the arrays atom_ptr_i and atom_ptr_j are defined in read_residuedata. The array indexing should be such that i_atoms[i] and atom_ptr_i[i] correspond to the same atom. */ var_var_ENERGY var_var_energy_calc(milli_pdbATOM *i_atoms, milli_pdbATOM *j_atoms, ATOMRESPARAM **atom_ptr_i, ATOMRESPARAM **atom_ptr_j, bool *do_these_hbond_flag, bool *i_hbond_status, bool *j_hbond_status) { double r2; /* squared distance between two atoms */ double E, EE; double xmax, ymax, zmax, xmin, ymin, zmin; var_var_ENERGY energy; double r, f, alpha; double dx, dy, dz; int i, j, disulf_flag; int i_hbond_index, j_hbond_index, donor_index_wrt_hydrogen; bool atompair_hbond_flag; extern int COULOMB_FLAG, GB_FLAG; VDW_ENERGY vdwE; extern double **CHARGE_PRODUCT, VDW_CUTOFF; *do_these_hbond_flag=0; E=0; energy.E_vdw = 0; energy.E_coulomb = 0; energy.E_pol = 0; VDW_ENERGY_0(vdwE); i_hbond_index=0; j_hbond_index=0; i=1; while(atom_ptr_i[i]!=NULL) { if(atom_ptr_i[i]->atom_ptr->atom_class!=ENDFLAG) { xmax=i_atoms[i].coord.x+FORCEFIELD_DISTANCE_CUTOFF; xmin=i_atoms[i].coord.x-FORCEFIELD_DISTANCE_CUTOFF; ymax=i_atoms[i].coord.y+FORCEFIELD_DISTANCE_CUTOFF; ymin=i_atoms[i].coord.y-FORCEFIELD_DISTANCE_CUTOFF; zmax=i_atoms[i].coord.z+FORCEFIELD_DISTANCE_CUTOFF; zmin=i_atoms[i].coord.z-FORCEFIELD_DISTANCE_CUTOFF; j=1; j_hbond_index=0; while(atom_ptr_j[j]!=NULL) { if(atom_ptr_j[j]->atom_ptr->atom_class!=ENDFLAG) { /* only measure distances that are within a cube of side 2*FORCEFIELD_DISTANCE_CUTOFF centered at j_atoms->coord */ if(j_atoms[j].coord.x<=xmax) if(j_atoms[j].coord.x>=xmin) if(j_atoms[j].coord.y<=ymax) if(j_atoms[j].coord.y>=ymin) if(j_atoms[j].coord.z<=zmax) if(j_atoms[j].coord.z>=zmin) { dx = j_atoms[j].coord.x - i_atoms[i].coord.x; dy = j_atoms[j].coord.y - i_atoms[i].coord.y; dz = j_atoms[j].coord.z - i_atoms[i].coord.z; r2 = dx*dx + dy*dy + dz*dz; // atoms on top of each other, yet both are real atoms // may be a problem with ligands if(r2<=TINY) vdwE.repulsive += VDW_CUTOFF; if(r2>TINY && r2<=FORCEFIELD_DISTANCE_CUTOFF_SQRD) { disulf_flag=0; /* don't score disulfide-bonded sulfurs */ if(r2<=MAX_DISULFIDE_BOND_LENGTH_SQRD) { if(strcmp(atom_ptr_i[i]->atom_ptr->atomtype,"SH")==0) if(strcmp(atom_ptr_j[j]->atom_ptr->atomtype,"SH")==0) disulf_flag=1; } if(disulf_flag==0) { r=sqrt(r2); if(atom_ptr_i[i]->charge!=0) if(atom_ptr_j[j]->charge!=0) /* if both are charged */ { if(GB_FLAG==1) { alpha = (i_atoms[i].born_radius)*(j_atoms[j].born_radius); f = f_of_r2_alpha(r2, alpha); EE = CHARGE_PRODUCT[atom_ptr_i[i]->charge_class][atom_ptr_j[j]->charge_class]/f; if(EE > POL_ATTRACTION_CAP) EE = POL_ATTRACTION_CAP; energy.E_pol += BORN_CONST(f)*EE; } if(COULOMB_FLAG==1) { energy.E_coulomb += coulomb_energy_q2_per_angstrom(r,atom_ptr_i[i],atom_ptr_j[j]); } } add_atompair_vdw_energy(&vdwE, atom_ptr_i[i], atom_ptr_j[j], r2, r); energy.E_hbond += rotamerpair_hbond_energy(&i_atoms[i], &j_atoms[j], &atom_ptr_i[i], &atom_ptr_j[j], r, r2, &atompair_hbond_flag, &donor_index_wrt_hydrogen); if(atompair_hbond_flag==1) // hbonding atompair { *do_these_hbond_flag=1; i_hbond_status[i_hbond_index]=1; j_hbond_status[j_hbond_index]=1; // mark donor if(atom_ptr_i[i]->donorflag[0] == 'H') i_hbond_status[i_hbond_index + donor_index_wrt_hydrogen]=1; else j_hbond_status[j_hbond_index + donor_index_wrt_hydrogen]=1; } } } } } if(atom_ptr_j[j]->donorflag[0] != ' ') ++j_hbond_index; ++j; } } if(atom_ptr_i[i]->donorflag[0] != ' ') ++i_hbond_index; ++i; } energy.E_coulomb = ((double)COULOMB_FLAG)*COULOMB_CONST_OVER_INTERNAL_DIELECTRIC*energy.E_coulomb; energy.E_pol = COULOMB_CONST*energy.E_pol; energy.E_vdw = sum_VDW_ENERGY(vdwE); return energy; } #include "output_stuff.h" /* returns var_fix_ENERGY between sidechain atoms var_atoms and fixed_atoms, as well as born radii for var_atoms gb_flag 1 means calculate born radii and energies; 0 means do not */ var_fix_ENERGY var_fixed_energy_calc(mini_pdbATOM *fixed_atoms, mini_pdbATOM *var_atoms, int gb_flag) { double r2; /* squared distance between two atoms */ double E, EE; var_fix_ENERGY energy; double E_vdw_14, E_coulomb_14, E_torsion; mini_pdbATOM *varatom1, *fixatom1, *contactAtom, *twobondAtom; double xmax, ymax, zmax, xmin, ymin, zmin; double r, alpha, f; static COULOMBIC *first_coulombic = NULL; static COULOMBIC *coulombic = NULL; int flag2; int fix_ctr, var_ctr, step; double dx, dy, dz; ATOMRESPARAM *bkbn_atomresparam; extern ATOMRESPARAM PSEB; /* pseudoatom sidechain for Born radii precalculations */ extern double VAR_FIX_BORN_SCALE; extern int GB_FLAG, COULOMB_FLAG, TORSION_FLAG; VDW_ENERGY vdwE; extern double **CHARGE_PRODUCT; E=0; EE=0; contactAtom = NULL; twobondAtom = NULL; VAR_FIX_ENERGY_0(energy); VDW_ENERGY_0(vdwE); if(gb_flag==2) { if(first_coulombic != NULL) { coulombic = first_coulombic; freeCOULOMBIC(coulombic); coulombic=NULL; first_coulombic=NULL; } return(energy); } E_vdw_14=0; E_coulomb_14=0; E_torsion=0; if(gb_flag==1) if(GB_FLAG == 1) { if(first_coulombic == NULL) { coulombic = (COULOMBIC *)malloc(sizeof(COULOMBIC)); first_coulombic = coulombic; coulombic->next = NULL; coulombic->r2 = ENDFLAG; } coulombic = first_coulombic; coulombic->r2 = ENDFLAG; } ++var_atoms; /* advance pointer to var_atoms[1] */ varatom1 = var_atoms; while(var_atoms->atom_ptr!=NULL) { var_atoms->hbond_satisfied = 1; if(var_atoms->atom_ptr->donorflag[0] != ' ') var_atoms->hbond_satisfied = 0; ++var_atoms; } var_atoms = varatom1; ++fixed_atoms; /* advance pointer to fixed_atoms[1] */ fixatom1 = fixed_atoms; while(fixed_atoms->atom_ptr!=NULL) { fixed_atoms->hbond_satisfied = 1; if(fixed_atoms->atom_ptr->donorflag[0] != ' ') fixed_atoms->hbond_satisfied = 0; ++fixed_atoms; } fixed_atoms = fixatom1; fix_ctr=1; var_ctr=1; while(fixed_atoms->atom_ptr!=NULL) { if(fixed_atoms->atom_ptr->atom_ptr->atom_class!=ENDFLAG) /* no U */ { xmax=fixed_atoms->coord.x+FORCEFIELD_DISTANCE_CUTOFF; xmin=fixed_atoms->coord.x-FORCEFIELD_DISTANCE_CUTOFF; ymax=fixed_atoms->coord.y+FORCEFIELD_DISTANCE_CUTOFF; ymin=fixed_atoms->coord.y-FORCEFIELD_DISTANCE_CUTOFF; zmax=fixed_atoms->coord.z+FORCEFIELD_DISTANCE_CUTOFF; zmin=fixed_atoms->coord.z-FORCEFIELD_DISTANCE_CUTOFF; var_atoms = varatom1; var_ctr=1; while(var_atoms->atom_ptr!=NULL) { if(var_atoms->atom_ptr->atom_ptr->atom_class!=ENDFLAG) { /* only measure distances that are within a cube of side 2*FORCEFIELD_DISTANCE_CUTOFF centered at var_atoms->coord */ if(var_atoms->coord.x<=xmax) if(var_atoms->coord.x>=xmin) if(var_atoms->coord.y<=ymax) if(var_atoms->coord.y>=ymin) if(var_atoms->coord.z<=zmax) if(var_atoms->coord.z>=zmin) { dx = var_atoms->coord.x - fixed_atoms->coord.x; dy = var_atoms->coord.y - fixed_atoms->coord.y; dz = var_atoms->coord.z - fixed_atoms->coord.z; r2 = dx*dx + dy*dy + dz*dz; if(r2 >= TINY && r2<=FORCEFIELD_DISTANCE_CUTOFF_SQRD) /* FORCEFIELD_DISTANCE_CUTOFF cutoff */ { flag2=0; E=0; r=sqrt(r2); bkbn_atomresparam = fixed_atoms->atom_ptr; /* find atom connectivity */ if(var_atoms->seq_position==fixed_atoms->seq_position) /* intra-residue */ { if(fixed_atoms->atom_ptr->other_info<0) /* backbone-sidechain interactions only */ { if(strcmp(fixed_atoms->atom_ptr->atom_ptr->atomtype, "PSE")==0) flag2=4; else { /* >=1-4 nonbonded interactions w/ backbone */ if(strcmp(var_atoms->atom_ptr->atomname,"CB")==0) { if(strcmp(fixed_atoms->atom_ptr->atom_ptr->atomtype,"HN")==0 || strstr(fixed_atoms->atom_ptr->atomname,"O")!=0 ) { flag2=1; contactAtom = fixed_atoms; find_backbone_atom(contactAtom, "CA", 1); if(strcmp(fixed_atoms->atom_ptr->atom_ptr->atomtype,"HN")==0) { twobondAtom = fixed_atoms; find_backbone_atom(twobondAtom, "N", -1); } else { twobondAtom = fixed_atoms; find_backbone_atom(twobondAtom, "C", -1); } } } else /* HB or XG atoms NOT HG */ if((strstr(var_atoms->atom_ptr->atomname,"G")!=NULL && strstr(var_atoms->atom_ptr->atomname,"HG")==NULL) || strstr(var_atoms->atom_ptr->atomname,"HB")!=NULL) { if(strcmp(fixed_atoms->atom_ptr->atomname,"N")==0 || strcmp(fixed_atoms->atom_ptr->atomname,"C")==0 || strstr(fixed_atoms->atom_ptr->atomname,"HA")!=NULL) { flag2=1; twobondAtom = fixed_atoms; find_backbone_atom(twobondAtom, "CA", 1); contactAtom = var_atoms; find_sidechain_atom(contactAtom, 8, -1); /* CB */ } else if(strcmp(fixed_atoms->atom_ptr->atomname,"CA")!=0) flag2=2; } else { /* HG or XD atoms NOT HD */ if((strstr(var_atoms->atom_ptr->atomname,"D")!=NULL && strstr(var_atoms->atom_ptr->atomname,"HD")==NULL) || strstr(var_atoms->atom_ptr->atomname,"HG")!=NULL) { if(strcmp(fixed_atoms->atom_ptr->atomname,"CA")==0) { flag2=1; twobondAtom = var_atoms; /* twobond atom in this case is always CB */ find_sidechain_atom(twobondAtom, 8, -1); if(var_atoms->atom_ptr->ringDihedAtom!=fixed_atoms->atom_ptr->atomnumber) { if(var_atoms->atom_ptr->contactAtomNumberatom_ptr->atomnumber) step=-1; else /* find contact atom */ step=1; contactAtom = var_atoms; find_sidechain_atom(contactAtom, var_atoms->atom_ptr->contactAtomNumber, step); } else /* 1,4 interaction via ring closure */ { if(var_atoms->atom_ptr->ringContactAtomatom_ptr->atomnumber) step=-1; else step=1; contactAtom = var_atoms; find_sidechain_atom(contactAtom, var_atoms->atom_ptr->ringContactAtom, step); } } else flag2=2; } else { /* >1-4 interactions w/ backbone */ flag2=2; } } } } } else /* inter-residue */ { if(var_atoms->atom_ptr->atomnumber==8) /* CB */ { if(abs(fixed_atoms->seq_position - var_atoms->seq_position) > 1) flag2=2; else { if(fixed_atoms->seq_position == (var_atoms->seq_position+1)) { if(strcmp(fixed_atoms->atom_ptr->atomname, "N")==0) { flag2=1; twobondAtom = fixed_atoms; while(twobondAtom->seq_position!=var_atoms->seq_position) --twobondAtom; while(strcmp(twobondAtom->atom_ptr->atomname, "C")!=0) --twobondAtom; contactAtom = twobondAtom; while(strcmp(contactAtom->atom_ptr->atomname, "CA")!=0) --contactAtom; } else flag2=2; } else { if(fixed_atoms->seq_position == (var_atoms->seq_position-1)) { if(strcmp(fixed_atoms->atom_ptr->atomname, "C")==0) { flag2=1; twobondAtom = fixed_atoms; while(twobondAtom->seq_position!=var_atoms->seq_position) ++twobondAtom; while(strcmp(twobondAtom->atom_ptr->atomname, "N")!=0) ++twobondAtom; contactAtom = twobondAtom; while(strcmp(contactAtom->atom_ptr->atomname, "CA")!=0) ++contactAtom; } else flag2=2; } } } } else { flag2=2; } if(strcmp(fixed_atoms->atom_ptr->atom_ptr->atomtype, "PSE")==0) { bkbn_atomresparam = &PSEB; flag2=3; } } if(var_atoms->atom_ptr->charge!=0) /* if var_atom is charged */ { if(flag2<=3) /* Born factors only for >=1,4 interactions */ /* scaled Born factor for charged var_atoms */ if(fixed_atoms->atom_ptr->other_info<0) /* only use backbone atoms for calculating the born radii of var atoms */ { if(strcmp(fixed_atoms->atom_ptr->atom_ptr->atomtype, "HU")!=0) if(gb_flag==1) if(GB_FLAG==1) var_atoms->born_radius += VAR_FIX_BORN_SCALE*born_nonbond_calc(var_atoms->atom_ptr, bkbn_atomresparam, r2); } if(fixed_atoms->atom_ptr->charge!=0) /* if both are charged */ { if(gb_flag==1) if(GB_FLAG == 1) { coulombic->index_i = var_ctr; coulombic->index_j = fix_ctr; coulombic->r2 = r2; coulombic->flag = flag2; if(coulombic->next == NULL) { coulombic->next = (COULOMBIC *)malloc(sizeof(COULOMBIC)); coulombic = coulombic->next; coulombic->next = NULL; /* end of list marker unless next element replaces it */ } else coulombic = coulombic->next; coulombic->r2 = ENDFLAG; /* end of list */ } } } if(flag2<=2) { if(flag2==2) { add_atompair_vdw_energy(&vdwE,fixed_atoms->atom_ptr, var_atoms->atom_ptr, r2, r); energy.E_coulomb += coulomb_energy_q2_per_angstrom(r, fixed_atoms->atom_ptr, var_atoms->atom_ptr); /* also marks hbonding atoms; will identify in calling function */ energy.E_hbond += hbond_energy(fixed_atoms, var_atoms, r, r2); } if(flag2==1) /* 1_4 interaction */ { E_vdw_14 += vdw_energy_14(fixed_atoms->atom_ptr, var_atoms->atom_ptr, r2); E_coulomb_14 += coulomb_energy_q2_per_angstrom(r, fixed_atoms->atom_ptr, var_atoms->atom_ptr); E_torsion += torsion_energy(var_atoms, contactAtom, twobondAtom, fixed_atoms); } } } } } ++var_atoms; ++var_ctr; } } ++fixed_atoms; ++fix_ctr; } /* GB energy */ energy.E_born = 0; energy.E_pol = 0; if(gb_flag==1) if(GB_FLAG == 1) { var_atoms = varatom1; while(var_atoms->atom_ptr!=NULL) { if(var_atoms->atom_ptr->charge!=0) if(var_atoms->atom_ptr->atom_ptr->atom_class!=ENDFLAG) { if(var_atoms->born_radius > 0) var_atoms->born_radius = -1.0; var_atoms->born_radius = -HALF_COULOMB_CONST/var_atoms->born_radius; EE = var_atoms->atom_ptr->charge2/var_atoms->born_radius; energy.E_born += BORN_CONST(var_atoms->born_radius)*EE; } ++var_atoms; } var_atoms = varatom1-1; fixed_atoms = fixatom1-1; coulombic = first_coulombic; while(coulombic->r2!=ENDFLAG) { alpha = var_atoms[coulombic->index_i].born_radius*fixed_atoms[coulombic->index_j].born_radius; f = f_of_r2_alpha(coulombic->r2, alpha); EE = CHARGE_PRODUCT[var_atoms[coulombic->index_i].atom_ptr->charge_class][fixed_atoms[coulombic->index_j].atom_ptr->charge_class]/f; if(coulombic->flag==2) if(EE > POL_ATTRACTION_CAP) EE = POL_ATTRACTION_CAP; energy.E_pol += BORN_CONST(f)*EE; coulombic = coulombic->next; } energy.E_born = HALF_COULOMB_CONST*energy.E_born; energy.E_pol = COULOMB_CONST*energy.E_pol; } var_atoms = varatom1-1; fixed_atoms = fixatom1-1; energy.E_coulomb = ((double)COULOMB_FLAG)*COULOMB_CONST_OVER_INTERNAL_DIELECTRIC*energy.E_coulomb; energy.E_1_4 = ((double)TORSION_FLAG)*WEIGHT_1_4*(NONBOND_FACTOR_1_4*(E_vdw_14 + COULOMB_CONST*E_coulomb_14) + E_torsion); energy.E_vdw = sum_VDW_ENERGY(vdwE); energy.E_hbond = HBOND_FUNCTION_FLAG*energy.E_hbond; return energy; } /* minimalist version of var_fixed_energy_calc for rotamerlet local minimization; returns var_fix_vdw_ENERGY */ var_fix_vdw_ENERGY var_fixed_vdw_screen(mini_pdbATOM *fixed_atoms, CARTESIAN *var_atoms, ATOMRESPARAM **var_atom_ptr, int seq_position) { double r2; /* squared distance between two atoms */ double E; var_fix_vdw_ENERGY var_fix_vdw_energy; double xmax, ymax, zmax, xmin, ymin, zmin; int flag2; int fix_ctr, var_ctr; double dx, dy, dz; VDW_ENERGY vdwE; double r=0; E=0; VDW_ENERGY_0(vdwE); fix_ctr=1; while(fixed_atoms[fix_ctr].atom_ptr!=NULL) { if(fixed_atoms[fix_ctr].atom_ptr->atom_ptr->atom_class!=ENDFLAG) /* no U */ { xmax=fixed_atoms[fix_ctr].coord.x+FORCEFIELD_DISTANCE_CUTOFF; xmin=fixed_atoms[fix_ctr].coord.x-FORCEFIELD_DISTANCE_CUTOFF; ymax=fixed_atoms[fix_ctr].coord.y+FORCEFIELD_DISTANCE_CUTOFF; ymin=fixed_atoms[fix_ctr].coord.y-FORCEFIELD_DISTANCE_CUTOFF; zmax=fixed_atoms[fix_ctr].coord.z+FORCEFIELD_DISTANCE_CUTOFF; zmin=fixed_atoms[fix_ctr].coord.z-FORCEFIELD_DISTANCE_CUTOFF; var_ctr=1; while(var_atom_ptr[var_ctr]!=NULL) { if(var_atom_ptr[var_ctr]->atom_ptr->atom_class!=ENDFLAG) { /* only measure distances that are within a cube of side 2*FORCEFIELD_DISTANCE_CUTOFF centered at var_atoms->coord */ if(var_atoms[var_ctr].x<=xmax) if(var_atoms[var_ctr].x>=xmin) if(var_atoms[var_ctr].y<=ymax) if(var_atoms[var_ctr].y>=ymin) if(var_atoms[var_ctr].z<=zmax) if(var_atoms[var_ctr].z>=zmin) { dx = var_atoms[var_ctr].x - fixed_atoms[fix_ctr].coord.x; dy = var_atoms[var_ctr].y - fixed_atoms[fix_ctr].coord.y; dz = var_atoms[var_ctr].z - fixed_atoms[fix_ctr].coord.z; r2 = dx*dx + dy*dy + dz*dz; if(r2 >= TINY && r2<=FORCEFIELD_DISTANCE_CUTOFF_SQRD) /* FORCEFIELD_DISTANCE_CUTOFF cutoff */ { flag2=0; E=0; if(seq_position==fixed_atoms[fix_ctr].seq_position) /* intra-residue */ { if(fixed_atoms[fix_ctr].atom_ptr->other_info<0) /* backbone-sidechain interactions only */ { if(strcmp(fixed_atoms[fix_ctr].atom_ptr->atom_ptr->atomtype, "PSE")==0) flag2=4; else { /* >=1-4 nonbonded interactions w/ backbone */ if(strcmp(var_atom_ptr[var_ctr]->atomname,"CB")==0) { if(strcmp(fixed_atoms[fix_ctr].atom_ptr->atom_ptr->atomtype,"HN")==0 || strstr(fixed_atoms[fix_ctr].atom_ptr->atomname,"O")!=0 ) { flag2=1; } } else /* HB or XG atoms NOT HG */ if((strstr(var_atom_ptr[var_ctr]->atomname,"G")!=NULL && strstr(var_atom_ptr[var_ctr]->atomname,"HG")==NULL) || strstr(var_atom_ptr[var_ctr]->atomname,"HB")!=NULL) { if(strcmp(fixed_atoms[fix_ctr].atom_ptr->atomname,"N")==0 || strcmp(fixed_atoms[fix_ctr].atom_ptr->atomname,"C")==0 || strstr(fixed_atoms[fix_ctr].atom_ptr->atomname,"HA")!=NULL) { flag2=1; } else if(strcmp(fixed_atoms[fix_ctr].atom_ptr->atomname,"CA")!=0) flag2=2; } else { /* HG or XD atoms NOT HD */ if((strstr(var_atom_ptr[var_ctr]->atomname,"D")!=NULL && strstr(var_atom_ptr[var_ctr]->atomname,"HD")==NULL) || strstr(var_atom_ptr[var_ctr]->atomname,"HG")!=NULL) { if(strcmp(fixed_atoms[fix_ctr].atom_ptr->atomname,"CA")==0) { flag2=1; } else flag2=2; } else { /* >1-4 interactions w/ backbone */ flag2=2; } } } } } else /* inter-residue */ { if(var_atom_ptr[var_ctr]->atomnumber==8) /* CB */ { if(abs(fixed_atoms[fix_ctr].seq_position - seq_position) > 1) flag2=2; else { if(fixed_atoms[fix_ctr].seq_position == (seq_position+1)) { if(strcmp(fixed_atoms[fix_ctr].atom_ptr->atomname, "N")==0) { flag2=1; } else flag2=2; } else { if(fixed_atoms[fix_ctr].seq_position == (seq_position-1)) { if(strcmp(fixed_atoms[fix_ctr].atom_ptr->atomname, "C")==0) { flag2=1; } else flag2=2; } } } } else flag2=2; if(strcmp(fixed_atoms[fix_ctr].atom_ptr->atom_ptr->atomtype, "PSE")==0) { flag2=3; } } if(flag2!=0) { if(flag2==2) { r=sqrt(r2); add_atompair_vdw_energy(&vdwE, fixed_atoms[fix_ctr].atom_ptr, var_atom_ptr[var_ctr], r2, r); } } } } } ++var_ctr; } } ++fix_ctr; } var_fix_vdw_energy.E_vdw = sum_VDW_ENERGY(vdwE); return var_fix_vdw_energy; } /* version of var_var_vdw_calc for rotamerlet local minimization; returns E_vdw */ double var_var_vdw_screen(milli_pdbATOM *i_atoms, milli_pdbATOM *j_atoms, ATOMRESPARAM **atom_ptr_i, ATOMRESPARAM **atom_ptr_j) { double r2; /* squared distance between two atoms */ double xmax, ymax, zmax, xmin, ymin, zmin; double dx, dy, dz; var_var_ENERGY var_var_energy; int i, j, flag; VDW_ENERGY vdwE; double r=0; VDW_ENERGY_0(vdwE); i=1; while(atom_ptr_i[i]!=NULL) { if(atom_ptr_i[i]->atom_ptr->atom_class!=ENDFLAG) { xmax=i_atoms[i].coord.x+FORCEFIELD_DISTANCE_CUTOFF; xmin=i_atoms[i].coord.x-FORCEFIELD_DISTANCE_CUTOFF; ymax=i_atoms[i].coord.y+FORCEFIELD_DISTANCE_CUTOFF; ymin=i_atoms[i].coord.y-FORCEFIELD_DISTANCE_CUTOFF; zmax=i_atoms[i].coord.z+FORCEFIELD_DISTANCE_CUTOFF; zmin=i_atoms[i].coord.z-FORCEFIELD_DISTANCE_CUTOFF; j=1; while(atom_ptr_j[j]!=NULL) { if(atom_ptr_j[j]->atom_ptr->atom_class!=ENDFLAG) { /* only measure distances that are within a cube of side 2*FORCEFIELD_DISTANCE_CUTOFF centered at j_atoms->coord */ if(j_atoms[j].coord.x<=xmax) if(j_atoms[j].coord.x>=xmin) if(j_atoms[j].coord.y<=ymax) if(j_atoms[j].coord.y>=ymin) if(j_atoms[j].coord.z<=zmax) if(j_atoms[j].coord.z>=zmin) { dx = j_atoms[j].coord.x - i_atoms[i].coord.x; dy = j_atoms[j].coord.y - i_atoms[i].coord.y; dz = j_atoms[j].coord.z - i_atoms[i].coord.z; r2 = dx*dx + dy*dy + dz*dz; if(r2 >= TINY && r2<=FORCEFIELD_DISTANCE_CUTOFF_SQRD) { flag=1; if(r2<=MAX_DISULFIDE_BOND_LENGTH_SQRD) /* don't score disulfide-bonded sulfurs */ { if(strcmp(atom_ptr_i[i]->atom_ptr->atomtype,"SH")==0) if(strcmp(atom_ptr_j[j]->atom_ptr->atomtype,"SH")==0) flag=0; } if(flag==1) { r=sqrt(r2); add_atompair_vdw_energy(&vdwE, atom_ptr_i[i], atom_ptr_j[j], r2, r); } } } } ++j; } } ++i; } var_var_energy.E_vdw = sum_VDW_ENERGY(vdwE); return(var_var_energy.E_vdw); } /* version of var_var_vdw_screen for methyl local minimization; returns E_vdw */ double var_var_vdw_screen_micro(CARTESIAN *i_atoms, CARTESIAN *j_atoms, ATOMRESPARAM **atom_ptr_i, ATOMRESPARAM **atom_ptr_j) { double r2; /* squared distance between two atoms */ double xmax, ymax, zmax, xmin, ymin, zmin; double dx, dy, dz; int i, j, flag; var_var_ENERGY var_var_energy; double r=0; VDW_ENERGY vdwE; VDW_ENERGY_0(vdwE); i=1; while(atom_ptr_i[i]!=NULL) { if(atom_ptr_i[i]->atom_ptr->atom_class!=ENDFLAG) { xmax=i_atoms[i].x+FORCEFIELD_DISTANCE_CUTOFF; xmin=i_atoms[i].x-FORCEFIELD_DISTANCE_CUTOFF; ymax=i_atoms[i].y+FORCEFIELD_DISTANCE_CUTOFF; ymin=i_atoms[i].y-FORCEFIELD_DISTANCE_CUTOFF; zmax=i_atoms[i].z+FORCEFIELD_DISTANCE_CUTOFF; zmin=i_atoms[i].z-FORCEFIELD_DISTANCE_CUTOFF; j=1; while(atom_ptr_j[j]!=NULL) { if(atom_ptr_j[j]->atom_ptr->atom_class!=ENDFLAG) { /* only measure distances that are within a cube of side 2*FORCEFIELD_DISTANCE_CUTOFF centered at j_atoms->coord */ if(j_atoms[j].x<=xmax) if(j_atoms[j].x>=xmin) if(j_atoms[j].y<=ymax) if(j_atoms[j].y>=ymin) if(j_atoms[j].z<=zmax) if(j_atoms[j].z>=zmin) { dx = j_atoms[j].x - i_atoms[i].x; dy = j_atoms[j].y - i_atoms[i].y; dz = j_atoms[j].z - i_atoms[i].z; r2 = dx*dx + dy*dy + dz*dz; if(r2>=TINY && r2<=FORCEFIELD_DISTANCE_CUTOFF_SQRD) { flag=1; if(r2<=MAX_DISULFIDE_BOND_LENGTH_SQRD) /* don't score disulfide-bonded sulfurs */ { if(strcmp(atom_ptr_i[i]->atom_ptr->atomtype,"SH")==0) if(strcmp(atom_ptr_j[j]->atom_ptr->atomtype,"SH")==0) flag=0; } if(flag==1) { r=sqrt(r2); add_atompair_vdw_energy(&vdwE, atom_ptr_i[i], atom_ptr_j[j], r2, r); } } } } ++j; } } ++i; } var_var_energy.E_vdw = sum_VDW_ENERGY(vdwE); return(var_var_energy.E_vdw); } /* returns 1 if i and j have at least one atom for all rotamers within FORCEFIELD_DISTANCE_CUTOFF of each other */ int do_these_interact(milli_pdbATOM **i_atoms, milli_pdbATOM **j_atoms, int num_i, int num_j) { double r2; /* squared distance between two atoms */ double xmax, ymax, zmax, xmin, ymin, zmin; double dx, dy, dz; int i, j, ii, jj; ii=1; jj=0; while(i_atoms[ii]!=NULL) { for(i=num_i;i>=1;--i) /* tip of sidechain -> backbone; greater chance that tips interact */ { xmax=i_atoms[ii][i].coord.x+FORCEFIELD_DISTANCE_CUTOFF; xmin=i_atoms[ii][i].coord.x-FORCEFIELD_DISTANCE_CUTOFF; ymax=i_atoms[ii][i].coord.y+FORCEFIELD_DISTANCE_CUTOFF; ymin=i_atoms[ii][i].coord.y-FORCEFIELD_DISTANCE_CUTOFF; zmax=i_atoms[ii][i].coord.z+FORCEFIELD_DISTANCE_CUTOFF; zmin=i_atoms[ii][i].coord.z-FORCEFIELD_DISTANCE_CUTOFF; jj=1; while(j_atoms[jj]!=NULL) { for(j=num_j;j>=1;--j) { if(j_atoms[jj][j].coord.x<=xmax) if(j_atoms[jj][j].coord.x>=xmin) if(j_atoms[jj][j].coord.y<=ymax) if(j_atoms[jj][j].coord.y>=ymin) if(j_atoms[jj][j].coord.z<=zmax) if(j_atoms[jj][j].coord.z>=zmin) { dx = j_atoms[jj][j].coord.x - i_atoms[ii][i].coord.x; dy = j_atoms[jj][j].coord.y - i_atoms[ii][i].coord.y; dz = j_atoms[jj][j].coord.z - i_atoms[ii][i].coord.z; r2 = dx*dx + dy*dy + dz*dz; if(r2<=FORCEFIELD_DISTANCE_CUTOFF_SQRD) return(1); } } ++jj; } } ++ii; } return(0); } void trim_LOOKUP_ROTAMER_X_HBOND(LOOKUP_ROTAMER_X_HBOND **lookupRotX_hbond) { int h; LOOKUP_ROTAMER_X_HBOND *holding_lookupRotX_hbond; if((*lookupRotX_hbond) == NULL) return; if((*lookupRotX_hbond)->hbond_rot_pair == NULL) return; holding_lookupRotX_hbond = (*lookupRotX_hbond); (*lookupRotX_hbond) = (LOOKUP_ROTAMER_X_HBOND *)malloc(sizeof(LOOKUP_ROTAMER_X_HBOND)); (*lookupRotX_hbond)->num_hbond_rot_pairs = holding_lookupRotX_hbond->num_hbond_rot_pairs; (*lookupRotX_hbond)->hbond_rot_pair = (HBOND_ROTAMER_PAIR *)calloc((*lookupRotX_hbond)->num_hbond_rot_pairs,sizeof(HBOND_ROTAMER_PAIR)); for(h=0;h<(*lookupRotX_hbond)->num_hbond_rot_pairs;++h) (*lookupRotX_hbond)->hbond_rot_pair[h] = holding_lookupRotX_hbond->hbond_rot_pair[h]; free_memory(holding_lookupRotX_hbond); } /* This function calculates the energy between two sidechain rotamers represented by i_gene and j_gene, using structures precalculated and stored in the lookup table and linked appropriately */ double rotamer_pair_energy(GENE i_gene, GENE j_gene, LOOKUP_ROTAMER_X_HBOND **lookupRotX_hbond) { var_var_ENERGY var_var_ENERGy; double var_var_energy, pair_electr; int h, k, x, y; double vdw_energy, dummy_vdw_energy1; 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; bool do_these_hbond_flag; static bool *i_hbond_status=NULL; static bool *j_hbond_status=NULL; extern int **BASE_THREE_CONVERT; extern int LOCAL_MINIMIZATION_FLAG; /* , MINIMIZE_METHYL_FLAG; */ extern int GB_FLAG, COULOMB_FLAG; if(i_gene->seq_position > j_gene->seq_position) failure_report("ERROR pairwise_energy_calc.cpp: rotamer_pair_energy must have i_gene->seq_position < j_gene->seq_position","abort"); best_rotamerlet_i = ENDFLAG; best_rotamerlet_j = ENDFLAG; if(i_hbond_status == NULL) { i_hbond_status = (bool *)calloc(MAX_RES_SIZE+2,sizeof(bool)); j_hbond_status = (bool *)calloc(MAX_RES_SIZE+2,sizeof(bool)); } if(i_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms != 0) if(j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms != 0) { falsify_bool_array(i_hbond_status, 0, MAX_RES_SIZE); falsify_bool_array(j_hbond_status, 0, MAX_RES_SIZE); } do_these_hbond_flag=0; var_var_ENERGy = var_var_energy_calc(i_gene->lookupRot->sideAtoms, j_gene->lookupRot->sideAtoms, i_gene->choice_ptr->resparam_ptr->atom_ptr, j_gene->choice_ptr->resparam_ptr->atom_ptr, &do_these_hbond_flag,i_hbond_status, j_hbond_status); pair_electr = ((double)COULOMB_FLAG)*var_var_ENERGy.E_coulomb + ((double)GB_FLAG)*var_var_ENERGy.E_pol; if(do_these_hbond_flag == 1) { if((*lookupRotX_hbond)==NULL) { (*lookupRotX_hbond) = (LOOKUP_ROTAMER_X_HBOND *)malloc(sizeof(LOOKUP_ROTAMER_X_HBOND)); (*lookupRotX_hbond)->hbond_rot_pair = (HBOND_ROTAMER_PAIR *)calloc(j_gene->choice_ptr->resparam_ptr->rotamerlib_ptr->number_of_rotamers+1,sizeof(HBOND_ROTAMER_PAIR)); (*lookupRotX_hbond)->num_hbond_rot_pairs = 0; (*lookupRotX_hbond)->hbond_rot_pair[0].j_index = ENDFLAG; } ++(*lookupRotX_hbond)->num_hbond_rot_pairs; h=0; while((*lookupRotX_hbond)->hbond_rot_pair[h].j_index != ENDFLAG) ++h; (*lookupRotX_hbond)->hbond_rot_pair[h+1].j_index = ENDFLAG; (*lookupRotX_hbond)->hbond_rot_pair[h].j_index = (short int)(j_gene->lookupRot_index-1); (*lookupRotX_hbond)->hbond_rot_pair[h].i_hbond_status = (bool *)calloc(i_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms,sizeof(bool)); (*lookupRotX_hbond)->hbond_rot_pair[h].j_hbond_status = (bool *)calloc(j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms,sizeof(bool)); copy_bool_array((*lookupRotX_hbond)->hbond_rot_pair[h].i_hbond_status, 0, i_hbond_status, 0, i_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms); copy_bool_array((*lookupRotX_hbond)->hbond_rot_pair[h].j_hbond_status, 0, j_hbond_status, 0, j_gene->choice_ptr->resparam_ptr->num_hbonding_sidechain_atoms); } 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<=jRES->rotamerlib_ptr->num_rotamerlets;++k) { vdw_energy = var_var_vdw_screen_micro(i_gene->lookupRot->rotamerlet[h].coord, j_gene->lookupRot->rotamerlet[k].coord, iRES->atom_ptr, jRES->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; } } } // if(MINIMIZE_METHYL_FLAG == 1) if(iRES->num_methyl_groups!=0 || jRES->num_methyl_groups!=0) /* relevant only if at least one of them has a methyl */ { if(iRES->num_methyl_groups!=0 && jRES->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<=jRES->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], j_gene->lookupRot->rotamerlet[best_rotamerlet_j].methyl_group->MeH[h].coord[k], iRES->MeH_atom_ptr[x], jRES->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], j_gene->lookupRot->rotamerlet[best_rotamerlet_j].methyl_group->non_MeH, iRES->MeH_atom_ptr[x], jRES->non_MeH_atom_ptr); } } for(h=1;h<=jRES->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, j_gene->lookupRot->rotamerlet[best_rotamerlet_j].methyl_group->MeH[h].coord[k], iRES->non_MeH_atom_ptr, jRES->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, j_gene->lookupRot->rotamerlet[best_rotamerlet_j].methyl_group->non_MeH, iRES->non_MeH_atom_ptr, jRES->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<=jRES->num_methyl_confs; ++conf_j) { vdw_energy = 0; for(x=1;x<=iRES->num_methyl_groups;++x) { for(h=1;h<=jRES->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<=jRES->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 = &j_gene->lookupRot->rotamerlet[best_rotamerlet_j]; resparam_ptr_non_Me = jRES; } else { rotamerlet_Me = &j_gene->lookupRot->rotamerlet[best_rotamerlet_j]; resparam_ptr_Me = jRES; rotamerlet_non_Me = &i_gene->lookupRot->rotamerlet[best_rotamerlet_i]; resparam_ptr_non_Me = iRES; } /* nonMeH-other_sidechain energies */ nonMeH_nonMeH_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 = nonMeH_nonMeH_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; } } } } } var_var_energy = var_var_ENERGy.E_vdw + pair_electr; return(var_var_energy); } ENERGY sidechain_internal_energy(mini_pdbATOM *side_pdb, COULOMBIC *first_coulombic, COULOMBIC *coulombic) { int pdb_ctr, pdb_ctr2, flag; VDW_ENERGY vdwE; ENERGY energy; double r2, r; double E_vdw, E_coulomb, E_vdw_14, E_coulomb_14, E_torsion; mini_pdbATOM *firstAtom, *contactAtom, *twobondAtom, *dihedAtom; /* zero the energies */ ENERGY_0(energy); VDW_ENERGY_0(vdwE); E_vdw = 0; E_coulomb = 0; E_vdw_14 = 0; E_coulomb_14 = 0; E_torsion = 0; coulombic = first_coulombic; if(coulombic!=NULL) { while(coulombic->r2!=ENDFLAG) coulombic = coulombic->next; } firstAtom = NULL; contactAtom = NULL; twobondAtom = NULL; dihedAtom = NULL; pdb_ctr = 1; while(side_pdb[pdb_ctr].seq_position!=ENDFLAG) { if(side_pdb[pdb_ctr].atom_ptr->atom_ptr->atom_class != ENDFLAG) { pdb_ctr2 = pdb_ctr+1; while(side_pdb[pdb_ctr2].seq_position!=ENDFLAG) { if(side_pdb[pdb_ctr2].atom_ptr->atom_ptr->atom_class != ENDFLAG) { r2 = Distance_sqrd(side_pdb[pdb_ctr].coord, side_pdb[pdb_ctr2].coord); if(r2>=TINY && r2<=FORCEFIELD_DISTANCE_CUTOFF_SQRD) { r=sqrt(r2); flag=0; /* Determine connectivity */ flag = get_atom_connectivity(r2, &side_pdb[pdb_ctr], &side_pdb[pdb_ctr2], &firstAtom, &contactAtom, &twobondAtom, &dihedAtom); if(coulombic!=NULL) { /* store pair interactions in coulombic linked-list for born energy calcs; we don't have Born radii yet, since this is just an isolated rotamer. When the approximate born radii become available, go through this list and calculate the energies (in generate_lookup_table). */ if(flag<=3) if(side_pdb[pdb_ctr].atom_ptr->charge!=0) if(side_pdb[pdb_ctr2].atom_ptr->charge!=0) { coulombic->index_i = pdb_ctr; coulombic->index_j = pdb_ctr2; coulombic->r2 = r2; coulombic->flag = flag; if(coulombic->next == NULL) { coulombic->next = (COULOMBIC *)malloc(sizeof(COULOMBIC)); coulombic = coulombic->next; coulombic->next = NULL; } else coulombic = coulombic->next; coulombic->r2 = ENDFLAG; } } if(flag!=0) { if(r2>=TINY && r2<=FORCEFIELD_DISTANCE_CUTOFF_SQRD) { if(flag==2) /* vdw energy */ { add_atompair_vdw_energy(&vdwE, side_pdb[pdb_ctr].atom_ptr, side_pdb[pdb_ctr2].atom_ptr, r2, r); E_coulomb += coulomb_energy_q2_per_angstrom(r,side_pdb[pdb_ctr].atom_ptr,side_pdb[pdb_ctr2].atom_ptr); } if(flag==1) /* 1,4 energy */ { E_vdw_14 += vdw_energy_14(side_pdb[pdb_ctr].atom_ptr, side_pdb[pdb_ctr2].atom_ptr, r2);; E_coulomb_14 += coulomb_energy_q2_per_angstrom(r,side_pdb[pdb_ctr].atom_ptr,side_pdb[pdb_ctr2].atom_ptr); E_torsion += torsion_energy(firstAtom, contactAtom, twobondAtom, dihedAtom); } } } } } ++pdb_ctr2; } } ++pdb_ctr; } energy.E_vdw = sum_VDW_ENERGY(vdwE); energy.E_coulomb = COULOMB_CONST_OVER_INTERNAL_DIELECTRIC*E_coulomb; energy.E_1_4 = WEIGHT_1_4*(E_torsion + NONBOND_FACTOR_1_4*(E_vdw_14 + COULOMB_CONST*E_coulomb_14)); coulombic = first_coulombic; return(energy); } double rotamer_bkbn_powell_objective_function(double *delta_chi, POWELL *powell) { int i,j,k; double E_forcefield; static mini_pdbATOM *side_pdb=NULL; static SIDECHAIN *side=NULL; var_fix_ENERGY var_fix_energy; ENERGY internal_forcefield_energy; COULOMBIC *coulombic, *first_coulombic; static double *chi; extern int MAX_RES_SIZE, NUM_FUNCTION_CALLS; ++NUM_FUNCTION_CALLS; /* don't let the rotamer stray too far from the library value */ for(i=1;i<=powell->minimize_rotamer_datatype->gene->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi;++i) { if(fabs(delta_chi[i]) >= powell->minimize_rotamer_datatype->gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[powell->minimize_rotamer_datatype->gene->lookupRot_index].half_std_dev[i]) return(10000000*POW2(delta_chi[i])); } if(side==NULL) { side = (SIDECHAIN *)malloc(sizeof(SIDECHAIN)); side->atom=(micro_pdbATOM *)calloc(MAX_RES_SIZE,sizeof(micro_pdbATOM)); side_pdb = (mini_pdbATOM *)calloc(MAX_RES_SIZE,sizeof(mini_pdbATOM)); chi = (double *)calloc(10,sizeof(double)); } coulombic = NULL; first_coulombic = NULL; for(i=1;i<=powell->minimize_rotamer_datatype->gene->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi;++i) { chi[i] = powell->minimize_rotamer_datatype->gene->chi[i] + delta_chi[i]; } build_a_SideChain(side, &powell->minimize_rotamer_datatype->gene->varpos_ptr->bkbn, powell->minimize_rotamer_datatype->gene->choice_ptr->resparam_ptr, &powell->minimize_rotamer_datatype->gene->seq_position, chi); side_pdb[1] = *powell->minimize_rotamer_datatype->CB_atom; /* CB */ k=9; j=2; while(side->atom[k].atom_ptr!=NULL) { side_pdb[j].seq_position = side_pdb[1].seq_position; side_pdb[j].coord = side->atom[k].coord; side_pdb[j].atom_ptr = powell->minimize_rotamer_datatype->gene->choice_ptr->resparam_ptr->atom_ptr[j]; ++k; ++j; } side_pdb[j].seq_position = ENDFLAG; side_pdb[j].atom_ptr = NULL; internal_forcefield_energy = sidechain_internal_energy(side_pdb, first_coulombic, coulombic); if(powell->minimize_rotamer_datatype->fixed_atoms!=NULL) var_fix_energy = var_fixed_energy_calc(powell->minimize_rotamer_datatype->fixed_atoms, side_pdb, 0); else VAR_FIX_ENERGY_0(var_fix_energy); E_forcefield = internal_forcefield_energy.E_vdw + internal_forcefield_energy.E_coulomb + internal_forcefield_energy.E_1_4 + var_fix_energy.E_vdw + var_fix_energy.E_coulomb + var_fix_energy.E_1_4; return(E_forcefield); } void minimize_rotamer_on_backbone(PROTEIN *protein, GENE gene, mini_pdbATOM *fixed_atoms, mini_pdbATOM *CB_atom) { static POWELL powell; int original_logfile_flag,h, original_max_function_calls, num_iter; double weight_vdw, weight_1_4; time_t start_time; double min_vdw, vdw_attractive_factor, vdw_repulsive_factor, vdw_cutoff; double **vdw_linear_switchpoint; static mini_pdbATOM *side_pdb=NULL, *dummyside_pdb=NULL; static double *delta_chi=NULL; extern int MAX_RES_SIZE, LOGFILE_FLAG, NUM_FUNCTION_CALLS; extern double VDW_ATTRACTIVE_FACTOR, VDW_REPULSIVE_FACTOR, VDW_CUTOFF; extern double **VDW_LINEAR_SWITCHPOINT; /* CYD shouldn't move; ALA doesn't really move; ligand coords are defined as exact */ if(strcmp(gene->choice_ptr->resparam_ptr->residuetype,"ALA")==0 || strcmp(gene->choice_ptr->resparam_ptr->residuetype,"CYD")==0 || gene->choice_ptr->resparam_ptr->ligand_flag==1) return; if(fixed_atoms==NULL) if(gene->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi==1) return; gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[gene->lookupRot_index].chi = (double *)calloc(gene->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi+2,sizeof(double)); for(h=1;h<=gene->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi;++h) gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[gene->lookupRot_index].chi[h] = gene->chi[h]; gene->chi = gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[gene->lookupRot_index].chi; vdw_attractive_factor = VDW_ATTRACTIVE_FACTOR; vdw_repulsive_factor = VDW_REPULSIVE_FACTOR; weight_vdw = WEIGHT_VDW; weight_1_4 = WEIGHT_1_4; vdw_cutoff = VDW_CUTOFF; vdw_linear_switchpoint = VDW_LINEAR_SWITCHPOINT; VDW_ATTRACTIVE_FACTOR = 1.0; VDW_REPULSIVE_FACTOR = 1.0; WEIGHT_VDW = 1.0; WEIGHT_1_4 = 1.0; VDW_CUTOFF = 1e10; VDW_LINEAR_SWITCHPOINT = NULL; original_logfile_flag = LOGFILE_FLAG; LOGFILE_FLAG=0; original_max_function_calls = protein->parameters.max_function_calls; protein->parameters.max_function_calls = 100*gene->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi; num_iter = gene->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi + 1; NUM_FUNCTION_CALLS=0; powell.protein = protein; if(delta_chi==NULL) { side_pdb = (mini_pdbATOM *)calloc((MAX_RES_SIZE+3), sizeof(mini_pdbATOM)); dummyside_pdb = (mini_pdbATOM *)calloc((MAX_RES_SIZE+3), sizeof(mini_pdbATOM)); delta_chi = (double *)calloc(10,sizeof(double)); powell.minimize_rotamer_datatype = (MINIMIZE_ROTAMER_DATATYPE *)malloc(sizeof(MINIMIZE_ROTAMER_DATATYPE)); } for(h=1;h<=gene->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi;++h) delta_chi[h]=0; powell.minimize_rotamer_datatype->gene = gene; powell.minimize_rotamer_datatype->fixed_atoms = fixed_atoms; powell.minimize_rotamer_datatype->CB_atom = CB_atom; min_vdw = rotamer_bkbn_powell_objective_function(delta_chi, &powell); start_time = time(NULL); powell_minimization(delta_chi, gene->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi, num_iter, &min_vdw, &rotamer_bkbn_powell_objective_function, &powell, start_time); for(h=1;h<=gene->choice_ptr->resparam_ptr->rotamerlib_ptr->numChi;++h) gene->chi[h] += delta_chi[h]; GENE_to_mini_pdbATOM(&gene,dummyside_pdb); side_pdb[1] = *CB_atom; h=2; while(dummyside_pdb[h-1].seq_position!=ENDFLAG) { side_pdb[h] = dummyside_pdb[h-1]; ++h; } side_pdb[h].sasa=ENDFLAG; side_pdb[h].seq_position=ENDFLAG; side_pdb[h].atom_ptr=NULL; VDW_REPULSIVE_FACTOR = vdw_repulsive_factor; VDW_ATTRACTIVE_FACTOR = vdw_attractive_factor; WEIGHT_VDW = weight_vdw; WEIGHT_1_4 = weight_1_4; VDW_CUTOFF = vdw_cutoff; VDW_LINEAR_SWITCHPOINT = vdw_linear_switchpoint; gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[gene->lookupRot_index].coulombic = (COULOMBIC *)malloc(sizeof(COULOMBIC)); gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[gene->lookupRot_index].first_coulombic = gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[gene->lookupRot_index].coulombic; gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[gene->lookupRot_index].coulombic->next = NULL; gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[gene->lookupRot_index].coulombic->r2 = ENDFLAG; gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[gene->lookupRot_index].energy = sidechain_internal_energy(side_pdb, gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[gene->lookupRot_index].first_coulombic, gene->choice_ptr->resparam_ptr->rotamerlib_ptr->rotamer[gene->lookupRot_index].coulombic); LOGFILE_FLAG = original_logfile_flag; protein->parameters.max_function_calls = original_max_function_calls; return; }