/* EGAD: energy_functions.h 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 is the header for energy_functions.cpp Also includes several macros and inline functions required for energy evaluations. */ #ifndef energy_functions_header_flag #define energy_functions_header_flag #include #include #include #include #include #include #include #include "structure_types.h" #include "moremath.h" #include "SASA.h" #include "dihedral_cartesian.h" #include "pdbATOM_utilities.h" /* constants used by many functions dependent on this file */ extern double NONBOND_FACTOR_1_4; extern double GENERAL_ASP; extern double HYDROPHOBIC_ASP; extern double WATER_DIELECTRIC; extern double ONE_OVER_WATER_DIELECTRIC; extern double INTERNAL_DIELECTRIC; extern double ONE_OVER_INTERNAL_DIELECTRIC; extern double BORN_TAU; extern double COULOMB_CONST_OVER_INTERNAL_DIELECTRIC; extern double COUL_CONST_BORN_TAU; extern double HALF_COUL_CONST_BORN_TAU; extern double IONIC_STRENGTH; extern double KAPPA; extern double WEIGHT_VDW, WEIGHT_ELECTROSTATICS, WEIGHT_1_4; extern double VDW_ATTRACTIVE_FACTOR,VDW_REPULSIVE_FACTOR,OVERALL_ENERGY_SCALE; extern int HBOND_FUNCTION_FLAG,GB_FLAG,SASA_FLAG,TORSION_FLAG,COULOMB_FLAG; extern double FORCEFIELD_DISTANCE_CUTOFF, FORCEFIELD_DISTANCE_CUTOFF_SQRD; /* returns f; distance-dependent & Born radii-dependent dielectric model of Still (JACS 112:6127-6129 (1990)) */ #define f_of_r2_alpha(r2,a) (sqrt((r2)+(a))) /* Born const for solvation energies; KAPPA is the ionic-strength-dependent debye constant */ #define BORN_CONST(f) ( KAPPA == 0 ? BORN_TAU : (ONE_OVER_WATER_DIELECTRIC*(exp(KAPPA*(f))) - ONE_OVER_INTERNAL_DIELECTRIC) ) /* This function calculates all energies for the structure in pdb; returned as ENERGY. pdb must have atom_ptrs assigned and be terminated with atom_ptr=NULL, seq_position=ENDFLAG. electr_calc_flag=0: no born calculations electr_calc_flag=1: calculate born radii and electrostatics. electr_calc_flag=2: calculate electrostatics w/ born radii already in pdb[i].born_radius. electr_calc_flag=3: store born desolvation sums in pdb[i].born_radius, but don't calculate born/pol energy sasa_flag=1: calculate surface-area based energies. surface area components also stored in extern SASA_SUM SASA_SUM_TEMP; this can be accessed by the calling function */ ENERGY energy_calc(mini_pdbATOM *pdb, int electr_calc_flag, int input_sasa_flag); /* set ENERGY components to 0 */ #define ENERGY_0(energy) \ { \ energy.E_vdw=0; \ energy.TdS=0; \ energy.E_coulomb=0; \ energy.E_1_4=0; \ energy.E_born=0; \ energy.E_pol=0; \ energy.E_sasa=0; \ energy.E_hbond=0; \ energy.E_folded=0; \ energy.E_unfolded=0; \ energy.E_rss=0; \ energy.E_reference=0; \ energy.E_specificity=0; \ energy.E_solubility=0; \ energy.E_structure=0; \ energy.E_total=0; \ energy.pseudo_dG=0; \ energy.clashes=0; \ energy.num_hbond=0; \ energy.num_bb_hbond=0; \ energy.num_sb_hbond=0; \ energy.num_ss_hbond=0; \ energy.num_unsatisfied_hbond=0; \ energy.num_side_unsatisfied_hbond=0; \ energy.num_bkbn_unsatisfied_hbond=0; \ energy.dG_gap=0; \ } /* scales each energy component by scale */ inline ENERGY scale_energy(double scale, ENERGY c) { ENERGY a; a.E_vdw = scale*c.E_vdw; a.TdS = scale*c.TdS; a.E_coulomb = scale*c.E_coulomb; a.E_1_4 = scale*c.E_1_4; a.E_born = scale*c.E_born; a.E_pol = scale*c.E_pol; a.E_sasa = scale*c.E_sasa; a.E_hbond = scale*c.E_hbond; a.E_reference = scale*c.E_reference; a.E_rss = scale*c.E_rss; a.E_specificity = scale*c.E_specificity; a.E_solubility = scale*c.E_solubility; a.E_folded = scale*c.E_folded; a.E_unfolded = scale*c.E_unfolded; a.E_total = scale*c.E_total; a.pseudo_dG = scale*c.pseudo_dG; a.E_structure = scale*c.E_structure; a.dG_gap = scale*c.dG_gap; return(a); } /* set VAR_FIX_ENERGY components to 0 */ #define VAR_FIX_ENERGY_0(energy) \ { \ energy.E_vdw=0; \ energy.E_coulomb=0; \ energy.E_1_4=0; \ energy.E_born=0; \ energy.E_pol=0; \ energy.E_sasa=0; \ energy.E_hbond=0; \ } /* set VAR_VAR_ENERGY components to 0 */ #define VAR_VAR_ENERGY_0(energy) \ { \ energy.E_vdw=0; \ energy.E_coulomb=0; \ energy.E_pol=0; \ double E_hbond; \ } /* adds components of VAR_FIX_ENERGY b and VAR_FIX_ENERGY c into VAR_FIX_ENERGY a */ #define VAR_FIX_ENERGY_ADD(a, b, c) \ { \ a.E_vdw = b.E_vdw + c.E_vdw; \ a.E_coulomb = b.E_coulomb + c.E_coulomb; \ a.E_1_4 = b.E_1_4 + c.E_1_4; \ a.E_born = b.E_born + c.E_born; \ a.E_pol = b.E_pol + c.E_pol; \ a.E_sasa = b.E_sasa + c.E_sasa; \ a.E_hbond = b.E_hbond + c.E_hbond; \ } /* subtracts components of ENERGY c from ENERGY b into ENERGY a */ inline ENERGY ENERGY_SUBTRACT(ENERGY b, ENERGY c) { ENERGY a; a.E_vdw = b.E_vdw - c.E_vdw; a.TdS = b.TdS - c.TdS; a.E_coulomb = b.E_coulomb - c.E_coulomb; a.E_1_4 = b.E_1_4 - c.E_1_4; a.E_born = b.E_born - c.E_born; a.E_pol = b.E_pol - c.E_pol; a.E_sasa = b.E_sasa - c.E_sasa; a.E_hbond = b.E_hbond - c.E_hbond; a.E_reference = b.E_reference - c.E_reference; a.E_specificity = b.E_specificity - c.E_specificity; a.E_solubility = b.E_solubility - c.E_solubility; a.E_folded = b.E_folded - c.E_folded; a.E_unfolded = b.E_unfolded - c.E_unfolded; a.E_rss = b.E_rss - c.E_rss; a.E_structure = b.E_structure - c.E_structure; a.pseudo_dG = b.pseudo_dG - c.pseudo_dG; a.E_total = b.E_total - c.E_total; a.clashes = b.clashes - c.clashes ; a.num_hbond = b.num_hbond - c.num_hbond ; a.num_bb_hbond = b.num_bb_hbond - c.num_bb_hbond ; a.num_sb_hbond = b.num_sb_hbond - c.num_sb_hbond ; a.num_ss_hbond = b.num_ss_hbond - c.num_ss_hbond ; a.num_unsatisfied_hbond = b.num_unsatisfied_hbond - c.num_unsatisfied_hbond ; a.num_side_unsatisfied_hbond = b.num_side_unsatisfied_hbond - c.num_side_unsatisfied_hbond ; a.num_bkbn_unsatisfied_hbond = b.num_bkbn_unsatisfied_hbond - c.num_bkbn_unsatisfied_hbond ; a.dG_gap = b.dG_gap - c.dG_gap ; return(a); } /* adds components of ENERGY b and ENERGY c into ENERGY a */ inline ENERGY ENERGY_ADD(ENERGY b, ENERGY c) { ENERGY a; a.E_vdw = b.E_vdw + c.E_vdw; a.TdS = b.TdS + c.TdS; a.E_coulomb = b.E_coulomb + c.E_coulomb; a.E_1_4 = b.E_1_4 + c.E_1_4; a.E_born = b.E_born + c.E_born; a.E_pol = b.E_pol + c.E_pol; a.E_sasa = b.E_sasa + c.E_sasa; a.E_hbond = b.E_hbond + c.E_hbond; a.E_reference = b.E_reference + c.E_reference; a.E_specificity = b.E_specificity + c.E_specificity; a.E_solubility = b.E_solubility + c.E_solubility; a.E_rss = b.E_rss + c.E_rss; a.E_folded = b.E_folded + c.E_folded; a.E_unfolded = b.E_unfolded + c.E_unfolded; a.E_structure = b.E_structure + c.E_structure; a.pseudo_dG = b.pseudo_dG + c.pseudo_dG; a.E_total = b.E_total + c.E_total; a.clashes = b.clashes + c.clashes; a.num_hbond = b.num_hbond + c.num_hbond; a.num_bb_hbond = b.num_bb_hbond + c.num_bb_hbond ; a.num_sb_hbond = b.num_sb_hbond + c.num_sb_hbond ; a.num_ss_hbond = b.num_ss_hbond + c.num_ss_hbond ; a.num_unsatisfied_hbond = b.num_unsatisfied_hbond + c.num_unsatisfied_hbond; a.num_side_unsatisfied_hbond = b.num_side_unsatisfied_hbond + c.num_side_unsatisfied_hbond; a.num_bkbn_unsatisfied_hbond = b.num_bkbn_unsatisfied_hbond + c.num_bkbn_unsatisfied_hbond; a.dG_gap = b.dG_gap + c.dG_gap; return(a); } /* return sum of var_var_ENERGY components */ #define var_var_ENERGY_TOTAL(energy) (energy.E_vdw + ((double)COULOMB_FLAG)*energy.E_coulomb + ((double)GB_FLAG)*energy.E_pol + ((double)HBOND_FUNCTION_FLAG)*energy.E_hbond) /* return sum of var_fix_ENERGY components */ #define var_fix_ENERGY_TOTAL(energy) (energy.E_vdw + ((double)COULOMB_FLAG)*energy.E_coulomb + ((double)TORSION_FLAG)*energy.E_1_4 + ((double)GB_FLAG)*(energy.E_born + energy.E_pol) + ((double)SASA_FLAG)*energy.E_sasa + ((double)HBOND_FUNCTION_FLAG)*energy.E_hbond) /* return sum of ENERGY components */ #define ENERGY_TOTAL(energy) (energy.E_vdw + ((double)COULOMB_FLAG)*energy.E_coulomb + ((double)TORSION_FLAG)*energy.E_1_4 + ((double)GB_FLAG)*(energy.E_pol + energy.E_born) + ((double)SASA_FLAG)*energy.E_sasa + ((double)HBOND_FUNCTION_FLAG)*energy.E_hbond) /* return sum of ENERGY components, scaled by OVERALL_ENERGY_SCALE and component usage flags; used for final output */ #define ENERGY_TOTAL_SCALED(energy) (OVERALL_ENERGY_SCALE*(energy.E_vdw + ((double)COULOMB_FLAG)*energy.E_coulomb + ((double)TORSION_FLAG)*energy.E_1_4 + ((double)GB_FLAG)*(energy.E_pol + energy.E_born) + ((double)SASA_FLAG)*energy.E_sasa + ((double)HBOND_FUNCTION_FLAG)*energy.E_hbond)) /* van der waals energies */ /* set VDW_ENERGY to zero */ #define VDW_ENERGY_0(vdwE) \ { \ vdwE.repulsive=0; \ vdwE.attractive=0; \ vdwE.linear=0; \ vdwE.clashes=0; \ } /* returns lennard_jones vdw energy between atoms i and j seperated by distance_squared r2 */ inline double lennard_jones_energy(double r2, ATOMRESPARAM *i, ATOMRESPARAM *j) { double sigma6_over_r6; extern double **VDW_SIGMA6_SCALE, **VDW_WELLDEPTH; sigma6_over_r6 = VDW_SIGMA6_SCALE[i->atom_ptr->atom_class][j->atom_ptr->atom_class]/(r2*r2*r2); return(VDW_WELLDEPTH[i->atom_ptr->atom_class][j->atom_ptr->atom_class]*sigma6_over_r6*(sigma6_over_r6 - 1)); } /* return 2 for side-side ; 1 for side-bkbn ; 0 for bkbn-bkbn or side-internal */ inline int identify_side_side(ATOMRESPARAM *i, ATOMRESPARAM *j, int seqpos_i, int seqpos_j) { int i_cb_flag, j_cb_flag; if(seqpos_i == seqpos_j) /* within same position; must be side-bkbn, bkbn-bkbn, or side-internal */ { return(0); } if(seqpos_i != seqpos_j) /* not in same position */ { if(i->other_info>0 && j->other_info>0) /* side-side */ { return(2); } else /* CB-CB, CB-side, side-bkbn, or bkbn-bkbn */ { i_cb_flag=0; j_cb_flag=0; if(strcmp(i->atomname,"CB")==0) i_cb_flag=1; if(strcmp(j->atomname,"CB")==0) j_cb_flag=1; if(i->other_info<0 && j->other_info<0) /* CB-CB or bkbn-bkbn */ { if(i_cb_flag==1 && j_cb_flag==1) /* CB-CB -> side-side */ { return(2); } else /* bkbn-bkbn */ { return(0); } } else /* CB-side or side-bkbn */ { if(i_cb_flag==1 || j_cb_flag==1) /* CB-side -> side-side */ { return(2); } else /* side-bkbn */ { return(1); } } } } return(0); } /* calculate vdw energies, and add components to VDW_ENERGY vdwE */ inline void add_atompair_vdw_energy(VDW_ENERGY *vdwE, ATOMRESPARAM *i, ATOMRESPARAM *j, double r2, double r) { double E; extern double VDW_CUTOFF, VDW_CLASH_ENERGY; extern double **VDW_LINEAR_SLOPE, **VDW_LINEAR_INTERCEPT, **VDW_LINEAR_SWITCHPOINT; E=0; /* scale factors are multiplied by sum_VDW_ENERGY, which is called by the calling function */ /* this section is for the linearized repulsive region */ if(VDW_LINEAR_SWITCHPOINT!=NULL) { if(r > VDW_LINEAR_SWITCHPOINT[i->atom_ptr->atom_class][j->atom_ptr->atom_class]) /* LJ region */ { E = lennard_jones_energy(r2,i,j); if(E<0) vdwE->attractive += E; else vdwE->repulsive += E; } else /* linear region */ { E = VDW_LINEAR_SLOPE[i->atom_ptr->atom_class][j->atom_ptr->atom_class]*r + VDW_LINEAR_INTERCEPT[i->atom_ptr->atom_class][j->atom_ptr->atom_class]; vdwE->linear += E; } } else /* for conventional LJ */ { E = lennard_jones_energy(r2,i,j); /* the energy is really big; truncate to prevent overflow errors */ if(E>VDW_CUTOFF) E=VDW_CUTOFF; if(E<0) vdwE->attractive += E; else vdwE->repulsive += E; } if(E>VDW_CLASH_ENERGY) ++vdwE->clashes; } /* returns the sum of the attractive and repulsive components of vdw, scaled appropriately */ inline double sum_VDW_ENERGY(VDW_ENERGY vdwE) { extern double WEIGHT_VDW, VDW_ATTRACTIVE_FACTOR, VDW_REPULSIVE_FACTOR; return(WEIGHT_VDW*(VDW_ATTRACTIVE_FACTOR*vdwE.attractive + VDW_REPULSIVE_FACTOR*(vdwE.repulsive + vdwE.linear))); } /* returns coulombic interaction energy between atoms i and j seperated by distance r Has NOT been multiplied by COULOMB_CONST; value is in units of e^2/A^2 (e=electron charge, A = angstrom). */ inline double coulomb_energy_q2_per_angstrom(double r, ATOMRESPARAM *i, ATOMRESPARAM *j) { double E; extern double **CHARGE_PRODUCT; E = CHARGE_PRODUCT[i->charge_class][j->charge_class]/r; /* prevent favorable electrostatic interactions between over-lapping atoms from artefactually overwhelming the lennard-jones */ if(E < COULOMBIC_ATTRACTION_CAP) E = COULOMBIC_ATTRACTION_CAP; return(E); } /* returns vdw energy for atoms seperated by 3 bonds; must be scaled by NONBOND_FACTOR_1_4 in the calling function */ inline double vdw_energy_14(ATOMRESPARAM *i, ATOMRESPARAM *j, double r2) { double sigma6_over_r6; extern double **VDW_SIGMA6_14, **VDW_WELLDEPTH_14; sigma6_over_r6 = VDW_SIGMA6_14[i->atom_ptr->atom_class][j->atom_ptr->atom_class]/(r2*r2*r2); return(VDW_WELLDEPTH_14[i->atom_ptr->atom_class][j->atom_ptr->atom_class]*sigma6_over_r6*(sigma6_over_r6 - 1)); } /* finds the target sidechain atom with atomnumber target_atomnumber; atom1 should be set to some atom within the residue of interest. step -/+ 1; best guess for whether the target atom is before/after atom1 in the mini_pdbATOM array. upon exit, atom1 = pointer to the atom of interest. See energy_functions.cpp: get_atom_connectivity for usage examples. */ #define find_sidechain_atom(atom1, target_atomnumber, step) \ { \ int halt; \ mini_pdbATOM *atom_pointer; \ \ atom_pointer = atom1; \ \ halt=0; \ while(halt==0 && atom_pointer->seq_position==atom1->seq_position) \ { \ if(atom_pointer->atom_ptr->atomnumber==target_atomnumber) \ halt=1; \ else \ atom_pointer+=step; \ } \ if(halt==0) \ { \ atom_pointer = atom1; \ while(halt==0 && atom_pointer->seq_position==atom1->seq_position) \ { \ if(atom_pointer->atom_ptr->atomnumber==target_atomnumber) \ halt=1; \ else \ atom_pointer-=step; \ } \ } \ if(halt==0) \ { \ atom_pointer = atom1; \ while(halt==0 && atom_pointer->seq_position!=ENDFLAG) \ { \ if(atom_pointer->atom_ptr->atomnumber==target_atomnumber && \ atom_pointer->seq_position==atom1->seq_position) \ halt=1; \ else \ atom_pointer-=step; \ } \ } \ if(halt==0) \ { \ atom_pointer = atom1; \ while(halt==0 && atom_pointer->seq_position!=ENDFLAG) \ { \ if(atom_pointer->atom_ptr->atomnumber==target_atomnumber && \ atom_pointer->seq_position==atom1->seq_position) \ halt=1; \ else \ atom_pointer+=step; \ } \ } \ if(halt==0) \ { \ fprintf(stderr, "ERROR cannot find atom seq_position: %d atomnumber: %d\n", atom1->seq_position, target_atomnumber); \ abort(); \ } \ atom1 = atom_pointer; \ } /* finds the target backbone atom with atomname target_atomname; atom1 should be set to some atom within the residue of interest. step -/+ 1; best guess for whether the target atom is before/after atom1 in the mini_pdbATOM array. upon exit, atom1 = pointer to the atom of interest. See energy_functions.cpp: get_atom_connectivity for usage examples. */ #define find_backbone_atom(atom1, target_atomname, step) \ { \ int halt; \ mini_pdbATOM *atom_pointer; \ \ atom_pointer = atom1; \ \ halt=0; \ while(halt==0 && atom_pointer->seq_position==atom1->seq_position) \ { \ if(strcmp(atom_pointer->atom_ptr->atomname, target_atomname)==0) \ halt=1; \ else \ atom_pointer+=step; \ } \ if(halt==0) \ { \ atom_pointer = atom1; \ while(halt==0 && atom_pointer->seq_position==atom1->seq_position) \ { \ if(strcmp(atom_pointer->atom_ptr->atomname, target_atomname)==0) \ halt=1; \ else \ atom_pointer-=step; \ } \ } \ if(halt==0) \ { \ atom_pointer = atom1; \ while(halt==0 && atom_pointer->seq_position!=ENDFLAG) \ { \ if(strcmp(atom_pointer->atom_ptr->atomname, target_atomname)==0 \ && atom_pointer->seq_position==atom1->seq_position) \ halt=1; \ else \ atom_pointer-=step; \ } \ } \ if(halt==0) \ { \ atom_pointer = atom1; \ while(halt==0 && atom_pointer->seq_position!=ENDFLAG) \ { \ if(strcmp(atom_pointer->atom_ptr->atomname, target_atomname)==0 \ && atom_pointer->seq_position==atom1->seq_position) \ halt=1; \ else \ atom_pointer+=step; \ } \ } \ if(halt==0) \ { \ fprintf(stderr, "ERROR cannot find atom seq_position: %d name: %s\n", atom1->seq_position, target_atomname); \ exit(1); \ } \ atom1 = atom_pointer; \ } /* This function pre-calculates the components of desolvation due to atoms whose relative distances will not change during any of the calculations in EGAD (atoms <=1,3 and atoms on the same rigid ring). Called by read_forcefield.cpp: read_forcefield. GB implementation of Dominy and Brooks J Phys Chem B 103:3765-3773 (1999) */ void born_bonded(RESPARAM *inputresparam); /* This function precalculates the intrinsic energy for all rotamers. vdw energies between atoms of the same aromatic ring system are NOT calculated. 1,4 and bonded born factors are added and stored in rotamer->born_factor array. Electrostatics between charged atoms are measured and stored in COULOMBIC list. called by read_forcefield.cpp: read_forcefield and ligand_stuff.cpp: read_ligand_data.cpp */ void intrinsic_rotamer_energy(RESPARAM *resparam); /* Given the four atoms, return the dihedral-angle-dependent energy component of 1,4 interactions firstAtom \ contactAtom--------twobondAtom \ dihedAtom */ double torsion_energy(mini_pdbATOM *firstAtom, mini_pdbATOM *contactAtom, mini_pdbATOM *twobondAtom, mini_pdbATOM *dihedAtom); /* free a COULOMBIC linked-list; standard recursive linked-list freeing algorithm */ void freeCOULOMBIC(COULOMBIC *coulombic); /* calculates the desolvation components for a pair of non-bonded (>=1,4) atoms for Born radii calcs GB implementation of Dominy and Brooks J Phys Chem B 103:3765-3773 (1999) */ double born_nonbond_calc(ATOMRESPARAM *atom_i, ATOMRESPARAM *atom_j, double r2); /* determines connectivity between two atoms pdb_i and pdb_j in an array of mini_pdbATOM in the calling function returns 0 for <=1,3; 1 for 1,4; 2 for >1,4 ; for 1,4, returns non-NULL values for dihedral atoms */ int get_atom_connectivity(double r2, mini_pdbATOM *pdb_i, mini_pdbATOM *pdb_j, mini_pdbATOM **firstAtom, mini_pdbATOM **contactAtom, mini_pdbATOM **twobondAtom, mini_pdbATOM **dihedAtom); void build_hbond_burial_virtual_atoms(mini_pdbATOM *probe, mini_pdbATOM *pdb); double exposed_hbond_groups(mini_pdbATOM *pdb); #define find_donor_hydrogen(donor_hydrogen, donor_heavy) \ { \ donor_hydrogen = donor_heavy; \ \ while(donor_hydrogen->atom_ptr->donorflag[0]!='H' && \ donor_hydrogen->atom_ptr->contactAtomNumber != donor_heavy->atom_ptr->atomnumber && \ donor_hydrogen->seq_position == donor_heavy->seq_position) \ ++donor_hydrogen; \ \ if(donor_hydrogen->seq_position != donor_heavy->seq_position) \ { \ donor_hydrogen = donor_heavy; \ \ while(donor_hydrogen->atom_ptr->donorflag[0]!='H' && \ donor_hydrogen->atom_ptr->contactAtomNumber != donor_heavy->atom_ptr->atomnumber && \ donor_hydrogen->seq_position == donor_heavy->seq_position) \ --donor_hydrogen; \ } \ } #define find_donor_heavy(donor_heavy, donor_hydrogen) \ { \ donor_heavy = donor_hydrogen; \ \ while((donor_heavy->atom_ptr->donorflag[0]!='D' && \ donor_heavy->atom_ptr->donorflag[0]!='B') && \ donor_heavy->atom_ptr->atomnumber != donor_hydrogen->atom_ptr->contactAtomNumber && \ donor_heavy->seq_position == donor_hydrogen->seq_position) \ --donor_heavy; \ \ if(donor_hydrogen->seq_position != donor_heavy->seq_position) \ { \ donor_heavy = donor_hydrogen; \ \ while((donor_heavy->atom_ptr->donorflag[0]!='D' && \ donor_heavy->atom_ptr->donorflag[0]!='B') && \ donor_heavy->atom_ptr->atomnumber != donor_hydrogen->atom_ptr->contactAtomNumber && \ donor_heavy->seq_position == donor_hydrogen->seq_position) \ ++donor_heavy; \ } \ } inline bool can_these_hbond(ATOMRESPARAM *atom_1, ATOMRESPARAM *atom_2, double r) { /* atoms too far apart */ if(r > 2.5) return(0); /* atoms too close */ if(r <= 1.4) return(0); /* one of these is not an hbonding atom */ if(atom_1->donorflag[0] == ' ') return(0); if(atom_2->donorflag[0] == ' ') return(0); /* if one of the atoms is a donor, score the bond when the hydrogen comes along */ if(atom_1->donorflag[0] == 'D') return(0); if(atom_2->donorflag[0] == 'D') return(0); /* both are H */ if(atom_1->donorflag[0] == 'H' && atom_2->donorflag[0] == 'H') return(0); /* both are acceptor */ if(atom_1->donorflag[0] == 'A' && atom_2->donorflag[0] == 'A') return(0); if(atom_1->donorflag[0] == 'B' && atom_2->donorflag[0] == 'B') return(0); if(atom_1->donorflag[0] == 'A' && atom_2->donorflag[0] == 'B') return(0); if(atom_1->donorflag[0] == 'B' && atom_2->donorflag[0] == 'A') return(0); return(1); } typedef struct { ENERGY energy; VDW_ENERGY vdw_da; VDW_ENERGY vdw_daa; int n_hbond; double E_eps8; double E_eps1_da; double E_eps1_daa; double E_eps1_ha; double E_eps1_haa; double E_1012; } SUPER_ENERGY; inline double hbond_function(ATOMRESPARAM *donor_atomptr, ATOMRESPARAM *hydrogen_atomptr, ATOMRESPARAM *acceptor_atomptr, ATOMRESPARAM *acceptor_proximal_atomptr, double r_da2, double r_daa2, double r_ha2, double r_haa2, double cos_theta, double cos_phi) { double E_da, E_daa, E_ha, E_haa ; extern double **CHARGE_PRODUCT; extern double COULOMB_CONST_HBOND_INTERNAL_DIELECTRIC; extern int HBOND_FUNCTION_TYPE; extern double HBOND_WELLDEPTH, HBOND_REPULSIVE, HBOND_ATTRACTIVE; extern SUPER_ENERGY super_energy; if(HBOND_FUNCTION_TYPE == 3) // used for flipping asn/gln/his rotamers when loading the pdb template return(-1000); // Dreiding-type angle-dependent 12-10 potential if(HBOND_FUNCTION_TYPE == 2) return(HBOND_WELLDEPTH*(HBOND_REPULSIVE/r_da2 - HBOND_ATTRACTIVE)*cos_theta*cos_theta*cos_phi*cos_phi/pow(r_da2,5)); /* OPLS-AA with dielectric eps = HBOND_DIELECTRIC just for hydrogen, acceptor, donor, acceptor_proximal atoms E = (energy with eps=HBOND_DIELECTRIC) - (energy with eps=INTERNAL_DIELECTRIC) since these will be counted elsewhere for INTERNAL_DIELECTRIC; subtract to remove it Yes, it's a redundant calculation, but it makes life a lot simpler!! */ E_ha = CHARGE_PRODUCT[hydrogen_atomptr->charge_class][acceptor_atomptr->charge_class]/sqrt(r_ha2); if(E_ha < COULOMBIC_ATTRACTION_CAP) E_ha = COULOMBIC_ATTRACTION_CAP; return(COULOMB_CONST_HBOND_INTERNAL_DIELECTRIC*E_ha); E_da = CHARGE_PRODUCT[donor_atomptr->charge_class][acceptor_atomptr->charge_class]/sqrt(r_da2); if(E_da < COULOMBIC_ATTRACTION_CAP) E_da = COULOMBIC_ATTRACTION_CAP; E_daa = CHARGE_PRODUCT[donor_atomptr->charge_class][acceptor_proximal_atomptr->charge_class]/sqrt(r_daa2); if(E_daa < COULOMBIC_ATTRACTION_CAP) E_daa = COULOMBIC_ATTRACTION_CAP; E_haa = CHARGE_PRODUCT[hydrogen_atomptr->charge_class][acceptor_proximal_atomptr->charge_class]/sqrt(r_haa2); if(E_haa < COULOMBIC_ATTRACTION_CAP) E_haa = COULOMBIC_ATTRACTION_CAP; return( COULOMB_CONST*(E_da + E_daa)/2.0 + COULOMB_CONST*(E_ha + E_haa) - COULOMB_CONST*(E_da + E_daa + E_ha + E_haa)/8.0 ); add_atompair_vdw_energy(&super_energy.vdw_da, donor_atomptr, acceptor_atomptr, r_da2, sqrt(r_da2)); add_atompair_vdw_energy(&super_energy.vdw_daa, donor_atomptr, acceptor_proximal_atomptr, r_daa2, sqrt(r_daa2)); ++super_energy.n_hbond; super_energy.E_eps8 += COULOMB_CONST*(E_da + E_daa + E_ha + E_haa)/8.0; super_energy.E_eps1_da += COULOMB_CONST*E_da; super_energy.E_eps1_daa += COULOMB_CONST*E_daa; super_energy.E_eps1_ha += COULOMB_CONST*E_ha; super_energy.E_eps1_haa += COULOMB_CONST*E_haa; super_energy.E_1012 += HBOND_WELLDEPTH*(HBOND_REPULSIVE/r_da2 - HBOND_ATTRACTIVE)*cos_theta*cos_theta*cos_phi*cos_phi/pow(r_da2,5); /* if(DEBUG == 1) printf("%d\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", identify_side_side(hydrogen_atomptr, acceptor_atomptr, 1, 2), sqrt(r_ha2),arccos(cos_theta), arccos(cos_phi), COULOMB_CONST*(E_da + E_daa + E_ha + E_haa), super_energy.vdw_da.attractive + super_energy.vdw_daa.attractive, (super_energy.vdw_da.repulsive + super_energy.vdw_da.linear) + (super_energy.vdw_daa.repulsive + super_energy.vdw_daa.linear) ); */ } inline double rotamerpair_hbond_energy(milli_pdbATOM *atom_1, milli_pdbATOM *atom_2, ATOMRESPARAM **atomptr_1, ATOMRESPARAM **atomptr_2, double r, double r2, bool *do_these_hbond, int *donor_index_wrt_hydrogen) { double r_hd2, r_da2, r_haa2, r_aaa2, r_daa2, cos_theta, cos_phi, cos_omega; milli_pdbATOM *donor, *acceptor, *hydrogen, *acceptor_proximal; ATOMRESPARAM *donor_atomptr, *acceptor_atomptr, *hydrogen_atomptr, *acceptor_proximal_atomptr; *do_these_hbond=0; if(can_these_hbond((*atomptr_1), (*atomptr_2), r)==0) return(0); if((*atomptr_1)->donorflag[0] == 'H') /* atom_1 is hydrogen ; atom_2 is acceptor */ { hydrogen = atom_1; hydrogen_atomptr = *atomptr_1; acceptor = atom_2; acceptor_atomptr = *atomptr_2; } else /* atom_2 is hydrogen ; atom_1 is acceptor */ { hydrogen = atom_2; hydrogen_atomptr = *atomptr_2; acceptor = atom_1; acceptor_atomptr = *atomptr_1; } /* find donor atom; should be attached to the hydrogen */ donor_atomptr = hydrogen_atomptr; donor = hydrogen; (*donor_index_wrt_hydrogen) = 0; while(donor_atomptr->atomnumber != hydrogen_atomptr->contactAtomNumber && donor_atomptr->donorflag[0]!='D' && donor_atomptr->donorflag[0]!='B') { --donor_atomptr; --donor; --(*donor_index_wrt_hydrogen); if(donor_atomptr == NULL) { donor_atomptr = hydrogen_atomptr; donor = hydrogen; (*donor_index_wrt_hydrogen) = 0; while(donor_atomptr->atomnumber != hydrogen_atomptr->contactAtomNumber && donor_atomptr->donorflag[0]!='D' && donor_atomptr->donorflag[0]!='B') { ++donor_atomptr; ++donor; ++(*donor_index_wrt_hydrogen); } } } r_da2 = Distance_sqrd(donor->coord,acceptor->coord); if(r_da2 > 16.0) return(0); r_hd2 = POW2(hydrogen_atomptr->bondlength); cos_theta = (r2 + r_hd2 - r_da2)/(2*r*hydrogen_atomptr->bondlength); if(cos_theta>0) return(0); /* find acceptor_proximal atom attached to acceptor */ acceptor_proximal = acceptor; acceptor_proximal_atomptr = acceptor_atomptr; while(acceptor_proximal_atomptr->atomnumber!=acceptor_atomptr->contactAtomNumber) { --acceptor_proximal; --acceptor_proximal_atomptr; if(acceptor_proximal_atomptr == NULL) { acceptor_proximal = acceptor; acceptor_proximal_atomptr = acceptor_atomptr; while(acceptor_proximal_atomptr->atomnumber!=acceptor_atomptr->contactAtomNumber) { ++acceptor_proximal; ++acceptor_proximal_atomptr; } } } r_haa2 = Distance_sqrd(hydrogen->coord,acceptor_proximal->coord); r_aaa2 = POW2(acceptor_atomptr->bondlength); cos_phi = (r2 + r_aaa2 - r_haa2)/(2*r*acceptor_atomptr->bondlength); if(cos_phi>0) return(0); r_daa2 = Distance_sqrd(donor->coord,acceptor_proximal->coord); cos_omega = (r_da2 + r_aaa2 - r_daa2)/(2*sqrt(r_da2*r_aaa2)); if(cos_omega>0) return(0); *do_these_hbond=1; if(HBOND_FUNCTION_FLAG == 1) return(hbond_function(donor_atomptr, hydrogen_atomptr, acceptor_atomptr, acceptor_proximal_atomptr, r_da2, r_daa2, r2, r_haa2, cos_theta, cos_phi)); return(0); } /* criterion from MacDonald and Thornton JMB 238: 777 (1994) */ /* energy from Dreiding-type 10-12 potential */ inline double hbond_energy(mini_pdbATOM *atom_1, mini_pdbATOM *atom_2, double r, double r2) { mini_pdbATOM *donor, *acceptor, *hydrogen, *acceptor_proximal; double r_hd2, r_da2, r_haa2, r_aaa2, r_daa2, cos_theta, cos_phi, cos_omega; if(can_these_hbond(atom_1->atom_ptr, atom_2->atom_ptr, r)==0) return(0); if(atom_1->atom_ptr->donorflag[0] == 'H') /* atom_1 is hydrogen ; atom_2 is acceptor */ { hydrogen = atom_1; acceptor = atom_2; } else /* atom_2 is hydrogen ; atom_1 is acceptor */ { hydrogen = atom_2; acceptor = atom_1; } /* find donor atom; should be attached to the hydrogen */ donor = hydrogen; find_donor_heavy(donor, hydrogen); r_da2 = Distance_sqrd(donor->coord,acceptor->coord); if(r_da2 > 16.0) return(0); r_hd2 = POW2(hydrogen->atom_ptr->bondlength); cos_theta = (r2 + r_hd2 - r_da2)/(2*r*hydrogen->atom_ptr->bondlength); if(cos_theta>0) return(0); /* find acceptor_proximal atom attached to acceptor */ acceptor_proximal = acceptor; find_sidechain_atom(acceptor_proximal, acceptor->atom_ptr->contactAtomNumber, -1); r_haa2 = Distance_sqrd(hydrogen->coord,acceptor_proximal->coord); r_aaa2 = POW2(acceptor->atom_ptr->bondlength); cos_phi = (r2 + r_aaa2 - r_haa2)/(2*r*acceptor->atom_ptr->bondlength); if(cos_phi>0) return(0); r_daa2 = Distance_sqrd(donor->coord,acceptor_proximal->coord); cos_omega = (r_da2 + r_aaa2 - r_daa2)/(2*sqrt(r_da2*r_aaa2)); if(cos_omega>0) return(0); atom_1->hbond_satisfied = 1; atom_2->hbond_satisfied = 1; donor->hbond_satisfied = 1; return(hbond_function(donor->atom_ptr, hydrogen->atom_ptr, acceptor->atom_ptr, acceptor_proximal->atom_ptr, r_da2, r_daa2, r2, r_haa2, cos_theta, cos_phi)); } /* return the energy for atom index_of_atom_of_interest within pdb array; born radii and sasa must already be calculated; used for diagnostics */ ENERGY energy_for_atom_of_interest(mini_pdbATOM *pdb, int index_of_atom_of_interest, int electr_flag, int sasa_flag, int torsion_flag); inline double hbonding_status_dependent_specificity_energy(bool hbonding_status, ATOMRESPARAM *atom_ptr) { extern int HBOND_SPECIFICITY_FLAG; if(HBOND_SPECIFICITY_FLAG == 0) return(0); if(atom_ptr->donorflag[0] == ' ') /* not a hbonding atom */ return(0); if(hbonding_status == 1) /* hbonding is satisfied */ return(0); return(atom_ptr->E_hbond_specificity); /* hbonding is unsatisfied; there may be an alternative state */ } #endif