smooth correction in host

host/inc/processligand.h

@@ -147,12 +147,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
@@ -190,8 +200,20 @@ 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,

host/src/main.cpp

@@ -93,8 +93,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

host/src/processligand.cpp

@@ -1100,6 +1100,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,
@@ -1124,6 +1125,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
@@ -1793,6 +1821,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,
@@ -1914,3 +1944,182 @@ 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);
+}
+

host/src/processresult.cpp

@@ -60,7 +60,7 @@ void write_basic_info(FILE* fp, const Liganddata* ligand_ref, const Dockpars* my
 
 
 	fprintf(fp, "***********************************\n");
-	fprintf(fp, "**   OCLADOCK REPORT FILE   **\n");
+	fprintf(fp, "**      OCLADOCK REPORT FILE     **\n");
 	fprintf(fp, "***********************************\n\n\n");
 
 	//Writing out docking parameters
@@ -153,7 +153,7 @@ void write_basic_info_dlg(FILE* fp, const Liganddata* ligand_ref, const Dockpars
 
 
 	fprintf(fp, "**********************************************************\n");
-	fprintf(fp, "**       OCLADOCK AUTODOCKTOOLS-COMPATIBLE DLG FILE       **\n");
+	fprintf(fp, "**      OCLADOCK AUTODOCKTOOLS-COMPATIBLE DLG FILE      **\n");
 	fprintf(fp, "**********************************************************\n\n\n");
 
 	//Writing out docking parameters
@@ -302,8 +302,10 @@ void make_resfiles(float* final_population, float* energies, const Liganddata* l
 		//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