corrected grid forces for ligand with rotbond

.gitignore

@@ -4,6 +4,7 @@ initpop.txt
 *.xml
 input/albumin_dock/
 ocladock.wiki/
+bin/
 final_population_run*
 device/stringify_tmp
 KernelProgramBuildInfo.txt

Makefile

@@ -197,7 +197,11 @@ odock: check-env-all stringify $(SRC)
 
 
 # Example
-PDB     := 1ac8
+# 1ac8: for testing gradients of translation and rotation genes
+# 7cpa: for testing gradients of torsion genes (15 torsions) 
+# 3tmn: for testing gradients of torsion genes (1 torsion)
+
+PDB     := 3tmn
 NRUN    := 1
 POPSIZE := 10
 TESTNAME:= test

device/calcenergy.cl

@@ -22,6 +22,9 @@ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 */
 
 
+#define DEBUG_ENERGY
+
+
 #include "calcenergy_basic.h"
 
 // All related pragmas are in defines.h (accesible by host and device code)
@@ -53,6 +56,11 @@ void gpu_calc_energy(
 		    	__local float* calc_coords_z,
 		    	__local float* partial_energies,
 
+			#if defined (DEBUG_ENERGY)
+			__local float* partial_interE,
+			__local float* partial_intraE,
+			#endif
+
 	             __constant float* atom_charges_const,
                      __constant char*  atom_types_const,
                      __constant char*  intraE_contributors_const,
@@ -76,6 +84,12 @@ void gpu_calc_energy(
 {
 	partial_energies[get_local_id(0)] = 0.0f;
 
+	#if defined (DEBUG_ENERGY)
+	partial_interE[get_local_id(0)] = 0.0f;
+	partial_intraE[get_local_id(0)] = 0.0f;
+	#endif
+
+
 #if defined (IMPROVE_GRID)
 	// INTERMOLECULAR for-loop (intermediate results)
 	// It stores a product of two chars
@@ -276,6 +290,10 @@ void gpu_calc_energy(
 				                  || (y >= dockpars_gridsize_y-1)
 						  || (z >= dockpars_gridsize_z-1)){
 			partial_energies[get_local_id(0)] += 16777216.0f; //100000.0f;
+	
+			#if defined (DEBUG_ENERGY)
+			partial_interE[get_local_id(0)] += 16777216.0f;
+			#endif
 		}
 		else
 		{
@@ -380,6 +398,10 @@ void gpu_calc_energy(
 			// Calculating affinity energy
 			partial_energies[get_local_id(0)] += TRILININTERPOL(cube, weights);
 
+			#if defined (DEBUG_ENERGY)
+			partial_interE[get_local_id(0)] += TRILININTERPOL(cube, weights);
+			#endif
+
 			// Capturing electrostatic values
 			atom_typeid = dockpars_num_of_atypes;
 
@@ -444,6 +466,10 @@ void gpu_calc_energy(
 			// Calculating electrosatic energy
 			partial_energies[get_local_id(0)] += q * TRILININTERPOL(cube, weights);
 
+			#if defined (DEBUG_ENERGY)
+			partial_interE[get_local_id(0)] += q * TRILININTERPOL(cube, weights);
+			#endif
+
 			// Capturing desolvation values
 			atom_typeid = dockpars_num_of_atypes+1;
 
@@ -507,10 +533,36 @@ void gpu_calc_energy(
 
 			// Calculating desolvation energy
 			partial_energies[get_local_id(0)] += fabs(q) * TRILININTERPOL(cube, weights);
+
+			#if defined (DEBUG_ENERGY)
+			partial_interE[get_local_id(0)] += fabs(q) * TRILININTERPOL(cube, weights);
+			#endif
 		}
 
 	} // End atom_id for-loop (INTERMOLECULAR ENERGY)
 
+
+	#if defined (DEBUG_ENERGY)
+	barrier(CLK_LOCAL_MEM_FENCE);
+
+	if (get_local_id(0) == 0)
+	{
+		float energy_interE = partial_interE[0];
+
+		for (uint contributor_counter=1;
+		          contributor_counter<NUM_OF_THREADS_PER_BLOCK;
+		          contributor_counter++)
+		{
+			energy_interE += partial_interE[contributor_counter];
+		}
+		partial_interE[0] = energy_interE;
+		//printf("%-20s %-10.8f\n", "energy_interE: ", energy_interE);
+	}
+
+	barrier(CLK_LOCAL_MEM_FENCE);
+	#endif
+
+
 	// In paper: intermolecular and internal energy calculation
 	// are independent from each other, -> NO BARRIER NEEDED
   	// but require different operations,
@@ -559,33 +611,49 @@ void gpu_calc_energy(
 			// Calculating van der Waals / hydrogen bond term
 #if defined (NATIVE_PRECISION)
 			partial_energies[get_local_id(0)] += native_divide(VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance,12));
+
+
+			#if defined (DEBUG_ENERGY)
+			partial_intraE[get_local_id(0)] += native_divide(VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance,12));
+			#endif
+
 #elif defined (HALF_PRECISION)
-			pfor (uint rotation_counter = get_local_id(0);
+			for (uint rotation_counter = get_local_id(0);
 	          rotation_counter < dockpars_rotbondlist_length;
 	          rotation_counter+=NUM_OF_THREADS_PER_BLOCK)
-	{artial_energies[get_local_id(0)] += half_divide(VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],half_powr(atomic_distance,12));
+	{partial_energies[get_local_id(0)] += half_divide(VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],half_powr(atomic_distance,12));
 #else	// Full precision
 			partial_energies[get_local_id(0)] += VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid]/powr(atomic_distance,12);
 #endif
 
-			if (intraE_contributors_const[3*contributor_counter+2] == 1)	//H-bond
+			if (intraE_contributors_const[3*contributor_counter+2] == 1) {	//H-bond
 #if defined (NATIVE_PRECISION)
 				partial_energies[get_local_id(0)] -= native_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance,10));
+			#if defined (DEBUG_ENERGY)
+			partial_intraE[get_local_id(0)] -= native_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance,10));
+			#endif
+
+
 #elif defined (HALF_PRECISION)
 				partial_energies[get_local_id(0)] -= half_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],half_powr(atomic_distance,10));
 #else	// Full precision
 				partial_energies[get_local_id(0)] -= VWpars_BD_const[atom1_typeid*dockpars_num_of_atypes+atom2_typeid]/powr(atomic_distance,10);
 #endif
-
-			else	//van der Waals
+			}
+			else {	//van der Waals
 #if defined (NATIVE_PRECISION)
 				partial_energies[get_local_id(0)] -= native_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance,6));
+
+			#if defined (DEBUG_ENERGY)
+			partial_intraE[get_local_id(0)] -= native_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance,6));
+			#endif
+
 #elif defined (HALF_PRECISION)
 				partial_energies[get_local_id(0)] -= half_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],half_powr(atomic_distance,6));
 #else	// Full precision
 				partial_energies[get_local_id(0)] -= VWpars_BD_const[atom1_typeid*dockpars_num_of_atypes+atom2_typeid]/powr(atomic_distance,6);
 #endif
-
+			}
 			// Calculating electrostatic term
 
 /*
@@ -600,7 +668,24 @@ void gpu_calc_energy(
                                                              atomic_distance * (-8.5525f + half_divide(86.9525f,(1.0f + 7.7839f*half_exp(-0.3154f*atomic_distance))))
                                                              );
 #else	// Full precision
-				partial_energies[get_local_id(0)] += dockpars_coeff_elec*atom_charges_const[atom1_id]*atom_charges_const[atom2_id]/
+				partial_energies[get_local_	#if defined (DEBUG_ENERGY)
+	barrier(CLK_LOCAL_MEM_FENCE);
+
+	if (get_local_id(0) == 0)
+	{
+		float energy_interE = partial_interE[0];
+
+		for (uint contributor_counter=1;
+		          contributor_counter<NUM_OF_THREADS_PER_BLOCK;
+		          contributor_counter++)
+		{
+			energy_interE += partial_interE[contributor_counter];
+		}
+	}
+	printf("%-20s %-10.8f\n", "energy_interE: ", energy_interE);
+
+	barrier(CLK_LOCAL_MEM_FENCE);
+	#endifid(0)] += dockpars_coeff_elec*atom_charges_const[atom1_id]*atom_charges_const[atom2_id]/
 			                                       (atomic_distance*(-8.5525f + 86.9525f/(1.0f + 7.7839f*exp(-0.3154f*atomic_distance))));
 #endif
 */
@@ -611,6 +696,15 @@ void gpu_calc_energy(
                                                              dockpars_coeff_elec * atom_charges_const[atom1_id] * atom_charges_const[atom2_id],
                                                              atomic_distance * (DIEL_A + native_divide(DIEL_B,(1.0f + DIEL_K*native_exp(-DIEL_B_TIMES_H*atomic_distance))))
                                                              );
+
+			#if defined (DEBUG_ENERGY)
+			partial_intraE[get_local_id(0)] += native_divide (
+                                                             dockpars_coeff_elec * atom_charges_const[atom1_id] * atom_charges_const[atom2_id],
+                                                             atomic_distance * (DIEL_A + native_divide(DIEL_B,(1.0f + DIEL_K*native_exp(-DIEL_B_TIMES_H*atomic_distance))))
+                                                             );
+			#endif
+
+
 #elif defined (HALF_PRECISION)
         partial_energies[get_local_id(0)] += half_divide (
                                                              dockpars_coeff_elec * atom_charges_const[atom1_id] * atom_charges_const[atom2_id],
@@ -635,6 +729,15 @@ void gpu_calc_energy(
 					                       (dspars_S_const[atom2_typeid] +
 							       dockpars_qasp*fabs(atom_charges_const[atom2_id]))*dspars_V_const[atom1_typeid]) *
 					                       dockpars_coeff_desolv*native_exp(-atomic_distance*native_divide(atomic_distance,25.92f));
+
+			#if defined (DEBUG_ENERGY)
+			partial_intraE[get_local_id(0)] += ((dspars_S_const[atom1_typeid] +
+							       dockpars_qasp*fabs(atom_charges_const[atom1_id]))*dspars_V_const[atom2_typeid] +
+					                       (dspars_S_const[atom2_typeid] +
+							       dockpars_qasp*fabs(atom_charges_const[atom2_id]))*dspars_V_const[atom1_typeid]) *
+					                       dockpars_coeff_desolv*native_exp(-atomic_distance*native_divide(atomic_distance,25.92f));
+			#endif
+
 #elif defined (HALF_PRECISION)
 			partial_energies[get_local_id(0)] += ((dspars_S_const[atom1_typeid] +
 							       dockpars_qasp*fabs(atom_charges_const[atom1_id]))*dspars_V_const[atom2_typeid] +
@@ -666,6 +769,26 @@ void gpu_calc_energy(
 		}
 	}
 
+	#if defined (DEBUG_ENERGY)
+	barrier(CLK_LOCAL_MEM_FENCE);
+
+	if (get_local_id(0) == 0)
+	{
+		float energy_intraE = partial_intraE[0];
+
+		for (uint contributor_counter=1;
+		          contributor_counter<NUM_OF_THREADS_PER_BLOCK;
+		          contributor_counter++)
+		{
+			energy_intraE += partial_intraE[contributor_counter];
+		}
+		partial_intraE[0] = energy_intraE;
+		//printf("%-20s %-10.8f\n", "energy_intraE: ", energy_intraE);
+
+	}
+	barrier(CLK_LOCAL_MEM_FENCE);
+	#endif
+
 }
 
 #include "kernel1.cl"

device/kernel1.cl

@@ -67,6 +67,10 @@ gpu_calc_initpop(
 	__local float calc_coords_y[MAX_NUM_OF_ATOMS];
 	__local float calc_coords_z[MAX_NUM_OF_ATOMS];
 	__local float partial_energies[NUM_OF_THREADS_PER_BLOCK];
+	#if defined (DEBUG_ENERGY)
+	__local float partial_interE [NUM_OF_THREADS_PER_BLOCK];
+	__local float partial_intraE [NUM_OF_THREADS_PER_BLOCK];
+	#endif
 
 	// Copying genotype from global memory
 	event_t ev = async_work_group_copy(genotype,
@@ -105,6 +109,10 @@ gpu_calc_initpop(
 			calc_coords_y,
 			calc_coords_z,
 			partial_energies,
+			#if defined (DEBUG_ENERGY)
+			partial_interE,
+			partial_intraE,
+			#endif
 
 	                atom_charges_const,
 		        atom_types_const,

device/kernel3.cl

@@ -97,6 +97,11 @@ perform_LS(
 	__local float calc_coords_z[MAX_NUM_OF_ATOMS];
 	__local float partial_energies[NUM_OF_THREADS_PER_BLOCK];
 
+	#if defined (DEBUG_ENERGY)
+	__local float partial_interE [NUM_OF_THREADS_PER_BLOCK];
+	__local float partial_intraE [NUM_OF_THREADS_PER_BLOCK];
+	#endif
+
 	// Determining run ID and entity ID
 	// Initializing offspring genotype
 	if (get_local_id(0) == 0)
@@ -204,6 +209,10 @@ perform_LS(
 				calc_coords_y,
 				calc_coords_z,
 				partial_energies,
+				#if defined (DEBUG_ENERGY)
+				partial_interE,
+				partial_intraE,
+				#endif
 
 				atom_charges_const,
 				atom_types_const,
@@ -298,6 +307,10 @@ perform_LS(
 					calc_coords_y,
 					calc_coords_z,
 					partial_energies,
+					#if defined (DEBUG_ENERGY)
+					partial_interE,
+					partial_intraE,
+					#endif
 
 					atom_charges_const,
 					atom_types_const,

device/kernel4.cl

@@ -89,6 +89,10 @@ gpu_gen_and_eval_newpops(
 	__local float calc_coords_y[MAX_NUM_OF_ATOMS];
 	__local float calc_coords_z[MAX_NUM_OF_ATOMS];
 	__local float partial_energies[NUM_OF_THREADS_PER_BLOCK];
+	#if defined (DEBUG_ENERGY)
+	__local float partial_interE [NUM_OF_THREADS_PER_BLOCK];
+	__local float partial_intraE [NUM_OF_THREADS_PER_BLOCK];
+	#endif
 
 	// In this case this compute-unit is responsible for elitist selection
 	if ((get_group_id(0) % dockpars_pop_size) == 0) {
@@ -255,6 +259,10 @@ gpu_gen_and_eval_newpops(
 				calc_coords_y,
 				calc_coords_z,
 				partial_energies,
+				#if defined (DEBUG_ENERGY)
+				partial_interE,
+				partial_intraE,
+				#endif
 		
                                 atom_charges_const,
 	                        atom_types_const,

device/kernel_gradient.cl

@@ -45,6 +45,9 @@ gradient_minimizer(
     	     __constant float* rotbonds_moving_vectors_const,
              __constant float* rotbonds_unit_vectors_const,
              __constant float* ref_orientation_quats_const,
+	     __constant int*   rotbonds_const,
+	     __constant int*   rotbonds_atoms_const,
+	     __constant int*   num_rotating_atoms_per_rotbond_const,
     			// Specific gradient-minimizer args
     //  __global float* restrict dockpars_conformations_next,   // initial population
                                                                 // whose (some) entities (genotypes) are to be minimized
@@ -199,23 +202,67 @@ gradient_minimizer(
 	__local float calc_coords_z[MAX_NUM_OF_ATOMS];
 	__local float partial_energies[NUM_OF_THREADS_PER_BLOCK];
 
+	#if defined (DEBUG_ENERGY)
+	__local float partial_interE [NUM_OF_THREADS_PER_BLOCK];
+	__local float partial_intraE [NUM_OF_THREADS_PER_BLOCK];
+	#endif
 	// -----------------------------------------------------------------------------
 	// Perform gradient-descent iterations
+
+		float grid_center_x = 52.340;
+		float grid_center_y = 15.029;
+		float grid_center_z = -2.932;
+		float ligand_center_x = 52.22740741;
+		float ligand_center_y = 15.51751852;
+		float ligand_center_z = -2.40896296;
+
 	do {
-///*
+/*
 		// Specific input genotypes for a ligand with no rotatable bonds (1ac8).
 		// Translation genes must be expressed in grids in OCLADock (local_genotype [0|1|2]).
 		// However, for testing purposes, 
 		// we start using translation values in real space (Angstrom): {31.79575, 93.743875, 47.699875}
-		// Rotation genes are expresed in the Shoemake space: local_genotype [3|4|5]
+		// Rotation genes are exp5.96902604 / (2*math.pi)
+resed in the Shoemake space: local_genotype [3|4|5]
 		// xyz_gene_gridspace = gridcenter_gridspace + (input_gene_realspace - gridcenter_realspace)/gridsize
+
+		// 1ac8				
 		local_genotype[0] = 30 + (31.79575  - 31.924) / 0.375;
 		local_genotype[1] = 30 + (93.743875 - 93.444) / 0.375;
 		local_genotype[2] = 30 + (47.699875 - 47.924) / 0.375;
 		local_genotype[3] = 0.1f;
 		local_genotype[4] = 0.5f;
 		local_genotype[5] = 0.9f;
+*/		
+
+///*
+		// 3tmn
+
+		local_genotype[0] = 30 + (ligand_center_x - grid_center_x) / 0.375;
+		local_genotype[1] = 30 + (ligand_center_y - grid_center_y) / 0.375;
+		local_genotype[2] = 30 + (ligand_center_z - grid_center_z) / 0.375;
+		local_genotype[3] = 0.02000000f;
+		local_genotype[4] = 1.44513262f / PI_TIMES_2;
+		local_genotype[5] = 5.96902604f / PI_TIMES_2;
+		local_genotype[6] = 0.0f;
+		local_genotype[7] = 0.0f;
+		local_genotype[8] = 0.0f;
+		local_genotype[9] = 0.0f;
+		local_genotype[10] = 0.0f;
+		local_genotype[11] = 0.0f;
+		local_genotype[12] = 0.0f;
+		local_genotype[13] = 0.0f;
+		local_genotype[14] = 0.0f;
+		local_genotype[15] = 0.0f;
+		local_genotype[16] = 0.0f;
+		local_genotype[17] = 0.0f;
+		local_genotype[18] = 0.0f;
+		local_genotype[19] = 0.0f;
+		local_genotype[20] = 0.0f;
+
+
 //*/
+
 		// Calculating gradient
 		// =============================================================
 		gpu_calc_gradient(
@@ -255,7 +302,10 @@ gradient_minimizer(
 				ref_coords_z_const,
 				rotbonds_moving_vectors_const,
 				rotbonds_unit_vectors_const,
-				ref_orientation_quats_const
+				ref_orientation_quats_const,
+				rotbonds_const,
+				rotbonds_atoms_const,
+				num_rotating_atoms_per_rotbond_const
 			 	// Gradient-related arguments
 			 	// Calculate gradients (forces) for intermolecular energy
 			 	// Derived from autodockdev/maps.py
@@ -346,15 +396,41 @@ gradient_minimizer(
 
 	// -----------------------------------------------------------------------------
 
-///*
+/*
+	// 1ac8
 	local_genotype[0] = 30 + (31.79575  - 31.924) / 0.375;
 	local_genotype[1] = 30 + (93.743875 - 93.444) / 0.375;
 	local_genotype[2] = 30 + (47.699875 - 47.924) / 0.375;
 	local_genotype[3] = 0.1f;
 	local_genotype[4] = 0.5f;
 	local_genotype[5] = 0.9f;
-//*/
+*/
 
+///*
+	// 7cpa
+
+	local_genotype[0] = 30 + (ligand_center_x - grid_center_x) / 0.375;
+	local_genotype[1] = 30 + (ligand_center_y - grid_center_y) / 0.375;
+	local_genotype[2] = 30 + (ligand_center_z - grid_center_z) / 0.375;
+	local_genotype[3] = 0.02000000f;
+	local_genotype[4] = 1.44513262f / PI_TIMES_2;
+	local_genotype[5] = 5.96902604f / PI_TIMES_2;
+		local_genotype[6] = 0.0f;
+		local_genotype[7] = 0.0f;
+		local_genotype[8] = 0.0f;
+		local_genotype[9] = 0.0f;
+		local_genotype[10] = 0.0f;
+		local_genotype[11] = 0.0f;
+		local_genotype[12] = 0.0f;
+		local_genotype[13] = 0.0f;
+		local_genotype[14] = 0.0f;
+		local_genotype[15] = 0.0f;
+		local_genotype[16] = 0.0f;
+		local_genotype[17] = 0.0f;
+		local_genotype[18] = 0.0f;
+		local_genotype[19] = 0.0f;
+		local_genotype[20] = 0.0f;
+//*/
 
   	// Calculating energy
 	// =============================================================
@@ -381,6 +457,10 @@ gradient_minimizer(
 			calc_coords_y,
 			calc_coords_z,
 			partial_energies,
+			#if defined (DEBUG_ENERGY)
+			partial_interE,
+			partial_intraE,
+			#endif
 
 	                atom_charges_const,
 		        atom_types_const,
@@ -401,6 +481,14 @@ gradient_minimizer(
 
 	//barrier(CLK_LOCAL_MEM_FENCE);
 
+	#if defined (DEBUG_ENERGY)
+	if (get_local_id(0) == 0) {
+		printf("%-20s %-10.8f\n", "GRIDE: ", partial_interE[0]);
+		printf("%-20s %-10.8f\n", "INTRAE: ", partial_intraE[0]);
+		printf("\n");
+	}
+	#endif
+
 	#if defined (DEBUG_MINIMIZER)
 	if (get_local_id(0) == 0) {
 		printf("%-40s %f\n", "AFTER- GRADIENT - local_energy: ", local_energy);

host/inc/calcenergy.h

@@ -100,6 +100,18 @@ typedef struct
        float rotbonds_moving_vectors_const[3*MAX_NUM_OF_ROTBONDS];
        float rotbonds_unit_vectors_const  [3*MAX_NUM_OF_ROTBONDS];
        float ref_orientation_quats_const  [4*MAX_NUM_OF_RUNS];
+	// Added for calculating torsion-related gradients.
+	// Passing list of rotbond-atoms ids to the GPU.
+	// Contains the same information as processligand.h/Liganddata->rotbonds 
+       int   rotbonds [2*MAX_NUM_OF_ROTBONDS];	
+
+	// Contains the same information as processligand.h/Liganddata->atom_rotbonds
+	// "atom_rotbonds": array that contains the rotatable bonds - atoms assignment.
+	// If the element atom_rotbonds[atom index][rotatable bond index] is equal to 1,
+	// it means,that the atom must be rotated if the bond rotates. A 0 means the opposite.
+       int  rotbonds_atoms [MAX_NUM_OF_ATOMS * MAX_NUM_OF_ROTBONDS];
+
+       int  num_rotating_atoms_per_rotbond [MAX_NUM_OF_ROTBONDS];
 } kernelconstant;
 
 int prepare_const_fields_for_gpu(Liganddata* 	   myligand_reference,

host/src/calcenergy.cpp

@@ -110,6 +110,29 @@ int prepare_const_fields_for_gpu(Liganddata*     myligand_reference,
 	float rotbonds_moving_vectors[3*MAX_NUM_OF_ROTBONDS];
 	float rotbonds_unit_vectors[3*MAX_NUM_OF_ROTBONDS];
 	float ref_orientation_quats[4*MAX_NUM_OF_RUNS];
+
+	// Added for calculating torsion-related gradients.
+	// Passing list of rotbond-atoms ids to the GPU.
+	// Contains the same information as processligand.h/Liganddata->rotbonds		
+
+	// Each row corresponds to one rotatable bond of the ligand.
+	// The rotatable bond is described with the indexes of the
+	// two atoms which are connected to each other by the bond.
+	// The row index is equal to the index of the rotatable bond.
+	int   rotbonds [2*MAX_NUM_OF_ROTBONDS];
+
+	// Contains the same information as processligand.h/Liganddata->atom_rotbonds
+	// "atom_rotbonds": array that contains the rotatable bonds - atoms assignment.
+	// If the element atom_rotbonds[atom index][rotatable bond index] is equal to 1,
+	// it means,that the atom must be rotated if the bond rotates. A 0 means the opposite.
+
+	// "rotbonds_atoms"
+	int  rotbonds_atoms [MAX_NUM_OF_ATOMS*MAX_NUM_OF_ROTBONDS];
+
+	// Each entry corresponds to a rotbond_id
+	// The value of an entry indicates the number of atoms that rotate 
+	// along with that rotbond_id
+	int  num_rotating_atoms_per_rotbond [MAX_NUM_OF_ROTBONDS];
 	// ------------------------------
 
 	//charges and type id-s
@@ -267,6 +290,53 @@ int prepare_const_fields_for_gpu(Liganddata*     myligand_reference,
 		//printf("Precalculated quaternion for run %d: %f %f %f %f\n", i, ref_orientation_quats[4*i], ref_orientation_quats[4*i+1], ref_orientation_quats[4*i+2], ref_orientation_quats[4*i+3]);
 	}
 
+	// Added for calculating torsion-related gradients.
+	// Passing list of rotbond-atoms ids to the GPU.
+	// Contains the same information as processligand.h/Liganddata->rotbonds
+	for (i=0; i < myligand_reference->num_of_rotbonds; i++)
+	{
+		rotbonds [2*i]   = myligand_reference->rotbonds[i][0]; // id of first-atom
+		rotbonds [2*i+1] = myligand_reference->rotbonds[i][1]; // id of second atom
+	}
+
+	// Contains the same information as processligand.h/Liganddata->atom_rotbonds
+	// "atom_rotbonds": array that contains the rotatable bonds - atoms assignment.
+	// If the element atom_rotbonds[atom index][rotatable bond index] is equal to 1,
+	// it means,that the atom must be rotated if the bond rotates. A 0 means the opposite.
+
+
+	for (i=0; i<MAX_NUM_OF_ROTBONDS; i++)
+	{
+		num_rotating_atoms_per_rotbond [i] = 0;
+	}
+
+
+	int* intpoi;
+	//intpoi = rotbonds_atoms;
+
+	for (i=0; i < myligand_reference->num_of_rotbonds; i++)
+	{	
+		// Pointing to the mem area corresponding to a given rotbond
+		intpoi = rotbonds_atoms + MAX_NUM_OF_ATOMS*i;
+
+		for (j=0; j < myligand_reference->num_of_atoms; j++)
+		{
+			/*
+			rotbonds_atoms [MAX_NUM_OF_ATOMS*i+j] = myligand_reference->atom_rotbonds [j][i]; // 
+			*/
+			
+			// If an atom rotates with a rotbond, then
+			// add its atom-id to the entry corresponding to the rotbond-id.
+			// Also, count the number of atoms that rotate with a certain rotbond
+			if (myligand_reference->atom_rotbonds [j][i] == 1){
+				*intpoi = j;
+				intpoi++;
+				num_rotating_atoms_per_rotbond [i] ++;	
+			}	
+
+		}
+	}
+
 	int m;
 	for (m=0;m<MAX_NUM_OF_ATOMS;m++){ KerConst->atom_charges_const[m] = atom_charges[m]; }
 	for (m=0;m<MAX_NUM_OF_ATOMS;m++){ KerConst->atom_types_const[m]   = atom_types[m]; }
@@ -282,7 +352,12 @@ int prepare_const_fields_for_gpu(Liganddata*     myligand_reference,
 	for (m=0;m<3*MAX_NUM_OF_ROTBONDS;m++){ KerConst->rotbonds_moving_vectors_const[m]= rotbonds_moving_vectors[m]; }
 	for (m=0;m<3*MAX_NUM_OF_ROTBONDS;m++){ KerConst->rotbonds_unit_vectors_const[m]  = rotbonds_unit_vectors[m]; }
 	for (m=0;m<4*MAX_NUM_OF_RUNS;m++)    { KerConst->ref_orientation_quats_const[m]  = ref_orientation_quats[m]; }
-
+	// Added for calculating torsion-related gradients.
+	// Passing list of rotbond-atoms ids to the GPU.
+	// Contains the same information as processligand.h/Liganddata->rotbonds
+	for (m=0;m<2*MAX_NUM_OF_ROTBONDS;m++) { KerConst->rotbonds[m] = rotbonds[m]; }
+	for (m=0;m<MAX_NUM_OF_ATOMS*MAX_NUM_OF_ROTBONDS;m++) { KerConst->rotbonds_atoms[m] = rotbonds_atoms[m]; }
+	for (m=0;m<MAX_NUM_OF_ROTBONDS;m++) { KerConst->num_rotating_atoms_per_rotbond[m] = num_rotating_atoms_per_rotbond[m]; }
 	return 0;
 }