Commit aec9fbd7 authored by Leonardo Solis's avatar Leonardo Solis

smooth correction in host

parent 6219813b
......@@ -127,12 +127,22 @@ double calc_ddd_Mehler_Solmajer(double);
int is_H_bond(const char*, const char*);
#if 0
void print_ref_lig_energies_f(Liganddata,
Gridinfo,
const float*,
const float,
const float,
const float);
#endif
void print_ref_lig_energies_f(Liganddata,
const float,
Gridinfo,
const float*,
const float,
const float,
const float);
//////////////////////////////////
//float functions
......@@ -170,8 +180,19 @@ void calc_interE_peratom_f(const Gridinfo* mygrid,
float peratom_elec [MAX_NUM_OF_ATOMS],
int debug);
#if 0
float calc_intraE_f(const Liganddata* myligand,
float dcutoff,
char ignore_desolv,
const float scaled_AD4_coeff_elec,
const float AD4_coeff_desolv,
const float qasp,
int debug);
#endif
float calc_intraE_f(const Liganddata* myligand,
float dcutoff,
float smooth,
char ignore_desolv,
const float scaled_AD4_coeff_elec,
const float AD4_coeff_desolv,
......
......@@ -81,8 +81,20 @@ int main(int argc, char* argv[])
//------------------------------------------------------------
// Calculating energies of reference ligand if required
//------------------------------------------------------------
#if 0
if (mypars.reflig_en_reqired == 1)
print_ref_lig_energies_f(myligand_init, mygrid, floatgrids, mypars.coeffs.scaled_AD4_coeff_elec, mypars.coeffs.AD4_coeff_desolv, mypars.qasp);
#endif
if (mypars.reflig_en_reqired == 1) {
print_ref_lig_energies_f(myligand_init,
mypars.smooth,
mygrid,
floatgrids,
mypars.coeffs.scaled_AD4_coeff_elec,
mypars.coeffs.AD4_coeff_desolv,
mypars.qasp);
}
//------------------------------------------------------------
// Starting Docking
......
......@@ -1064,6 +1064,7 @@ int is_H_bond(const char* atype1, const char* atype2)
return 0;
}
#if 0
void print_ref_lig_energies_f(Liganddata myligand,
Gridinfo mygrid,
const float* fgrids,
......@@ -1088,6 +1089,33 @@ void print_ref_lig_energies_f(Liganddata myligand,
printf("Intermolecular energy of reference ligand: %lf\n",
calc_interE_f(&mygrid, &myligand, fgrids, 0, 0));
}
#endif
void print_ref_lig_energies_f(Liganddata myligand,
const float smooth,
Gridinfo mygrid,
const float* fgrids,
const float scaled_AD4_coeff_elec,
const float AD4_coeff_desolv,
const float qasp)
//The function calculates the energies of the ligand given in the first parameter,
//and prints them to the screen.
{
double temp_vec [3];
int i;
printf("Intramolecular energy of reference ligand: %lf\n",
calc_intraE_f(&myligand, 8, smooth, 0, scaled_AD4_coeff_elec, AD4_coeff_desolv, qasp, 0));
for (i=0; i<3; i++)
temp_vec [i] = -1*mygrid.origo_real_xyz [i];
move_ligand(&myligand, temp_vec);
scale_ligand(&myligand, (double) 1.0/mygrid.spacing);
printf("Intermolecular energy of reference ligand: %lf\n",
calc_interE_f(&mygrid, &myligand, fgrids, 0, 0));
}
//////////////////////////////////
//float functions
......@@ -1641,6 +1669,8 @@ void calc_interE_peratom_f(const Gridinfo* mygrid,
//return interE;
}
// OCLADock original host "calc_intraE_f" function
#if 0
float calc_intraE_f(const Liganddata* myligand,
float dcutoff,
char ignore_desolv,
......@@ -1762,4 +1792,181 @@ float calc_intraE_f(const Liganddata* myligand,
else
return (vW + el);
}
#endif
// OCLADock host "calc_intraE_f" function
// corrected after smoothing was added
float calc_intraE_f(const Liganddata* myligand,
float dcutoff,
float smooth,
char ignore_desolv,
const float scaled_AD4_coeff_elec,
const float AD4_coeff_desolv,
const float qasp,
int debug)
//The function calculates the intramolecular energy of the ligand given by the first parameter,
//and returns it as a double. The second parameter is the distance cutoff, if the third isn't 0,
//desolvation energy won't be included by the energy value, the fourth indicates if messages
//about partial results are required (if debug=1)
{
int atom_id1, atom_id2;
int type_id1, type_id2;
float dist;
int distance_id;
int smoothed_distance_id;
float vdW1, vdW2;
float s1, s2, v1, v2;
float vW, el, desolv;
//The following tables will contain the 1/r^6, 1/r^10, 1/r^12, W_el/(r*eps(r)) and W_des*exp(-r^2/(2sigma^2)) functions for
//distances 0.01:0.01:20.48 A
static char first_call = 1;
static float r_6_table [2048];
static float r_10_table [2048];
static float r_12_table [2048];
static float r_epsr_table [2048];
static float desolv_table [2048];
//The following arrays will contain the q1*q2 and qasp*abs(q) values for the ligand which is the input parameter when this
//function is called first time (it is supposed that the energy must always be calculated for this ligand only, that is, there
//is only one ligand during the run of the program...)
static float q1q2 [256][256];
static float qasp_mul_absq [256];
//when first call, calculating tables
if (first_call == 1)
{
calc_distdep_tables_f(r_6_table, r_10_table, r_12_table, r_epsr_table, desolv_table, scaled_AD4_coeff_elec, AD4_coeff_desolv);
calc_q_tables_f(myligand, qasp, q1q2, qasp_mul_absq);
first_call = 0;
}
vW = 0;
el = 0;
desolv = 0;
if (debug == 1)
printf("\n\n\nINTRAMOLECULAR ENERGY CALCULATION\n\n");
for (atom_id1=0; atom_id1<myligand->num_of_atoms-1; atom_id1++) //for each atom pair
for (atom_id2=atom_id1+1; atom_id2<myligand->num_of_atoms; atom_id2++)
{
if (myligand->intraE_contributors [atom_id1][atom_id2] == 1) //if they have to be included in intramolecular energy calculation
{ //the energy contribution has to be calculated
dist = distance(&(myligand->atom_idxyzq [atom_id1][1]), &(myligand->atom_idxyzq [atom_id2][1]));
#if 0
if (dist <= 1)
{
if (debug == 1)
printf("\n\nToo low distance (%lf) between atoms %d and %d\n", dist, atom_id1, atom_id2);
//return HIGHEST_ENERGY; //returning maximal value
dist = 1;
}
#endif
if (debug == 1)
{
printf("\n\nCalculating energy contribution of atoms %d and %d\n", atom_id1+1, atom_id2+1);
printf("Distance: %lf\n", dist);
}
// Adding smoothing
// Getting type ids
type_id1 = myligand->atom_idxyzq [atom_id1][0];
type_id2 = myligand->atom_idxyzq [atom_id2][0];
unsigned int atom1_type_vdw_hb = myligand->atom1_types_reqm [type_id1];
unsigned int atom2_type_vdw_hb = myligand->atom2_types_reqm [type_id2];
// Getting optimum pair distance (opt_distance) from reqm and reqm_hbond
// reqm: equilibrium internuclear separation
// (sum of the vdW radii of two like atoms (A)) in the case of vdW
// reqm_hbond: equilibrium internuclear separation
// (sum of the vdW radii of two like atoms (A)) in the case of hbond
float opt_distance;
if (is_H_bond(myligand->atom_types [type_id1], myligand->atom_types [type_id2]) != 0) //H-bond
{
opt_distance = myligand->reqm_hbond [atom1_type_vdw_hb] + myligand->reqm_hbond [atom2_type_vdw_hb];
}
else //normal van der Waals
{
opt_distance = 0.5f*(myligand->reqm [atom1_type_vdw_hb] + myligand->reqm [atom2_type_vdw_hb]);
}
// Getting smoothed distance
// smoothed_distance = function(dist, opt_distance)
float smoothed_distance;
float delta_distance = 0.5f*smooth;
if (dist <= (opt_distance - delta_distance)) {
smoothed_distance = dist + delta_distance;
}
else if (dist < (opt_distance + delta_distance)) {
smoothed_distance = opt_distance;
}
else { // else if (dist >= (opt_distance + delta_distance))
smoothed_distance = dist - delta_distance;
}
distance_id = (int) floor((100*dist) + 0.5) - 1; // +0.5: rounding, -1: r_xx_table [0] corresponds to r=0.01
if (distance_id < 0) {
distance_id = 0;
}
smoothed_distance_id = (int) floor((100*smoothed_distance) + 0.5) - 1; // +0.5: rounding, -1: r_xx_table [0] corresponds to r=0.01
if (smoothed_distance_id < 0) {
smoothed_distance_id = 0;
}
#if 0
if ((dist < dcutoff) && (dist < 20.48)) //but only if the distance is less than distance cutoff value and 20.48A (because of the tables)
#endif
if (dist < dcutoff) //but only if the distance is less than distance cutoff value
{
if (is_H_bond(myligand->atom_types [type_id1], myligand->atom_types [type_id2]) != 0) //H-bond
{
vdW1 = myligand->VWpars_C [type_id1][type_id2]*r_12_table [smoothed_distance_id];
vdW2 = myligand->VWpars_D [type_id1][type_id2]*r_10_table [smoothed_distance_id];
if (debug == 1)
printf("H-bond interaction = ");
}
else //normal van der Waals
{
vdW1 = myligand->VWpars_A [type_id1][type_id2]*r_12_table [smoothed_distance_id];
vdW2 = myligand->VWpars_B [type_id1][type_id2]*r_6_table [smoothed_distance_id];
if (debug == 1)
printf("van der Waals interaction = ");
}
}
s1 = (myligand->solpar [type_id1] + qasp_mul_absq [atom_id1]);
s2 = (myligand->solpar [type_id2] + qasp_mul_absq [atom_id2]);
v1 = myligand->volume [type_id1];
v2 = myligand->volume [type_id2];
if (debug == 1)
printf(" %lf, electrostatic = %lf, desolv = %lf\n", (vdW1 - vdW2), q1q2[atom_id1][atom_id2] * r_epsr_table [distance_id],
(s1*v2 + s2*v1) * desolv_table [distance_id]);
vW += vdW1 - vdW2;
el += q1q2[atom_id1][atom_id2] * r_epsr_table [distance_id];
desolv += (s1*v2 + s2*v1) * desolv_table [distance_id];
#if 0
}
#endif
}
}
if (debug == 1)
printf("\nFinal energies: van der Waals = %lf, electrostatic = %lf, desolvation = %lf, total = %lf\n\n", vW, el, desolv, vW + el + desolv);
if (ignore_desolv == 0)
return (vW + el + desolv);
else
return (vW + el);
}
......@@ -324,8 +324,10 @@ void make_resfiles(float* final_population,
//if (i==127)
// accurate_intraE [i] = calc_intraE(&temp_docked, 8, 0, mypars->coeffs.scaled_AD4_coeff_elec, mypars->coeffs.AD4_coeff_desolv, 1); //calculating the intramolecular energy
//else
#if 0
accurate_intraE[i] = calc_intraE_f(&temp_docked, 8, 0, mypars->coeffs.scaled_AD4_coeff_elec, mypars->coeffs.AD4_coeff_desolv, mypars->qasp, debug);
#endif
accurate_intraE[i] = calc_intraE_f(&temp_docked, 8, mypars->smooth, 0, mypars->coeffs.scaled_AD4_coeff_elec, mypars->coeffs.AD4_coeff_desolv, mypars->qasp, debug);
move_ligand(&temp_docked, mygrid->origo_real_xyz); //moving it according to grid location
entity_rmsds [i] = calc_rmsd(ligand_from_pdb, &temp_docked, mypars->handle_symmetry); //calculating rmds compared to original pdb file
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment