calcenergy.cl 29.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
/*

OCLADock, an OpenCL implementation of AutoDock 4.2 running a Lamarckian Genetic Algorithm
Copyright (C) 2017 TU Darmstadt, Embedded Systems and Applications Group, Germany. All rights reserved.

AutoDock is a Trade Mark of the Scripps Research Institute.

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

*/

lvs's avatar
lvs committed
24
//#define DEBUG_ENERGY_KERNEL
25

26
27
// No needed to be included as all kernel sources are stringified
#if 0
28
#include "calcenergy_basic.h"
29
#endif
Leonardo Solis's avatar
Leonardo Solis committed
30

lvs's avatar
lvs committed
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
typedef struct
{
       float atom_charges_const[MAX_NUM_OF_ATOMS];
       char  atom_types_const  [MAX_NUM_OF_ATOMS];
} kernelconstant_interintra;

typedef struct
{
       char  intraE_contributors_const[3*MAX_INTRAE_CONTRIBUTORS];
} kernelconstant_intracontrib;

typedef struct
{
       float reqm_const [ATYPE_NUM];
       float reqm_hbond_const [ATYPE_NUM];
       unsigned int  atom1_types_reqm_const [ATYPE_NUM];
       unsigned int  atom2_types_reqm_const [ATYPE_NUM];
       float VWpars_AC_const   [MAX_NUM_OF_ATYPES*MAX_NUM_OF_ATYPES];
       float VWpars_BD_const   [MAX_NUM_OF_ATYPES*MAX_NUM_OF_ATYPES];
       float dspars_S_const    [MAX_NUM_OF_ATYPES];
       float dspars_V_const    [MAX_NUM_OF_ATYPES];
} kernelconstant_intra;

typedef struct
{
       int   rotlist_const     [MAX_NUM_OF_ROTATIONS];
} kernelconstant_rotlist;

typedef struct
{
       float ref_coords_x_const[MAX_NUM_OF_ATOMS];
       float ref_coords_y_const[MAX_NUM_OF_ATOMS];
       float ref_coords_z_const[MAX_NUM_OF_ATOMS];
       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];
} kernelconstant_conform;

Leonardo Solis's avatar
Leonardo Solis committed
69
70
// All related pragmas are in defines.h (accesible by host and device code)

71
72
73
74
75
76
void gpu_calc_energy(	    
				int    dockpars_rotbondlist_length,
				char   dockpars_num_of_atoms,
			    	char   dockpars_gridsize_x,
			    	char   dockpars_gridsize_y,
			    	char   dockpars_gridsize_z,
77
78
79
								    		// g1 = gridsize_x
				uint   dockpars_gridsize_x_times_y, 		// g2 = gridsize_x * gridsize_y
				uint   dockpars_gridsize_x_times_y_times_z,	// g3 = gridsize_x * gridsize_y * gridsize_z
80
81
82
83
84
85
86
		 __global const float* restrict dockpars_fgrids, // This is too large to be allocated in __constant 
		            	char   dockpars_num_of_atypes,
		            	int    dockpars_num_of_intraE_contributors,
			    	float  dockpars_grid_spacing,
			    	float  dockpars_coeff_elec,
			    	float  dockpars_qasp,
			    	float  dockpars_coeff_desolv,
lvs's avatar
lvs committed
87
				float  dockpars_smooth,
Leonardo Solis's avatar
Leonardo Solis committed
88

89
90
91
92
                    // Some OpenCL compilers don't allow declaring 
		    // local variables within non-kernel functions.
		    // These local variables must be declared in a kernel, 
		    // and then passed to non-kernel functions.
93
94
95
96
97
98
99
100
101
		    	__local float* genotype,
		   	__local float* energy,
		    	__local int*   run_id,

		    	__local float* calc_coords_x,
		    	__local float* calc_coords_y,
		    	__local float* calc_coords_z,
		    	__local float* partial_energies,

lvs's avatar
lvs committed
102
			#if defined (DEBUG_ENERGY_KERNEL)
103
104
105
			__local float* partial_interE,
			__local float* partial_intraE,
			#endif
lvs's avatar
lvs committed
106
107
108
#if 0
 				bool   debug,
#endif
lvs's avatar
lvs committed
109
110
111
112
113
		   __constant     kernelconstant_interintra* 		kerconst_interintra,
		   __global const kernelconstant_intracontrib*  	kerconst_intracontrib,
		   __constant     kernelconstant_intra*			kerconst_intra,
		   __constant     kernelconstant_rotlist*   		kerconst_rotlist,
		   __constant     kernelconstant_conform*		kerconst_conform
Leonardo Solis's avatar
Leonardo Solis committed
114
115
116
117
118
119
120
121
122
)

//The GPU device function calculates the energy of the entity described by genotype, dockpars and the liganddata
//arrays in constant memory and returns it in the energy parameter. The parameter run_id has to be equal to the ID
//of the run whose population includes the current entity (which can be determined with blockIdx.x), since this
//determines which reference orientation should be used.
{
	partial_energies[get_local_id(0)] = 0.0f;

lvs's avatar
lvs committed
123
	#if defined (DEBUG_ENERGY_KERNEL)
124
125
126
127
	partial_interE[get_local_id(0)] = 0.0f;
	partial_intraE[get_local_id(0)] = 0.0f;
	#endif

128
129
130
131
132
133
134
135
136
137
138
	// Convert orientation genes from sex. to radians
	float phi         = genotype[3] * DEG_TO_RAD;
	float theta       = genotype[4] * DEG_TO_RAD;
	float genrotangle = genotype[5] * DEG_TO_RAD;

	float genrot_unitvec [3];
	float sin_angle = native_sin(theta);
	genrot_unitvec [0] = sin_angle*native_cos(phi);
	genrot_unitvec [1] = sin_angle*native_sin(phi);
	genrot_unitvec [2] = native_cos(theta);

139
	uchar g1 = dockpars_gridsize_x;
lvs's avatar
lvs committed
140
141
	uint  g2 = dockpars_gridsize_x_times_y;
  	uint  g3 = dockpars_gridsize_x_times_y_times_z;
142

Leonardo Solis's avatar
Leonardo Solis committed
143
	// ================================================
Leonardo Solis's avatar
Leonardo Solis committed
144
	// CALCULATING ATOMIC POSITIONS AFTER ROTATIONS
Leonardo Solis's avatar
Leonardo Solis committed
145
	// ================================================
146
147
148
	for (uint rotation_counter = get_local_id(0);
	          rotation_counter < dockpars_rotbondlist_length;
	          rotation_counter+=NUM_OF_THREADS_PER_BLOCK)
Leonardo Solis's avatar
Leonardo Solis committed
149
	{
lvs's avatar
lvs committed
150
		int rotation_list_element = kerconst_rotlist->rotlist_const[rotation_counter];
Leonardo Solis's avatar
Leonardo Solis committed
151

152
		if ((rotation_list_element & RLIST_DUMMY_MASK) == 0)	// If not dummy rotation
Leonardo Solis's avatar
Leonardo Solis committed
153
		{
154
155
156
157
			uint atom_id = rotation_list_element & RLIST_ATOMID_MASK;

			// Capturing atom coordinates
			float atom_to_rotate[3];
Leonardo Solis's avatar
Leonardo Solis committed
158

159
			if ((rotation_list_element & RLIST_FIRSTROT_MASK) != 0)	// If first rotation of this atom
Leonardo Solis's avatar
Leonardo Solis committed
160
			{
lvs's avatar
lvs committed
161
162
163
				atom_to_rotate[0] = kerconst_conform->ref_coords_x_const[atom_id];
				atom_to_rotate[1] = kerconst_conform->ref_coords_y_const[atom_id];
				atom_to_rotate[2] = kerconst_conform->ref_coords_z_const[atom_id];
Leonardo Solis's avatar
Leonardo Solis committed
164
165
166
167
168
169
170
171
			}
			else
			{
				atom_to_rotate[0] = calc_coords_x[atom_id];
				atom_to_rotate[1] = calc_coords_y[atom_id];
				atom_to_rotate[2] = calc_coords_z[atom_id];
			}

172
			// Capturing rotation vectors and angle
173
			float rotation_unitvec[3];
174
			float rotation_movingvec[3];
175
			float rotation_angle;
176
177
178
179
180

			float quatrot_left_x, quatrot_left_y, quatrot_left_z, quatrot_left_q;
			float quatrot_temp_x, quatrot_temp_y, quatrot_temp_z, quatrot_temp_q;

			if ((rotation_list_element & RLIST_GENROT_MASK) != 0)	// If general rotation
Leonardo Solis's avatar
Leonardo Solis committed
181
			{
182
183
184
				rotation_unitvec[0] = genrot_unitvec[0];
				rotation_unitvec[1] = genrot_unitvec[1];
				rotation_unitvec[2] = genrot_unitvec[2];
185

Leonardo Solis's avatar
Leonardo Solis committed
186
187
188
				rotation_movingvec[0] = genotype[0];
				rotation_movingvec[1] = genotype[1];
				rotation_movingvec[2] = genotype[2];
189
190

				rotation_angle = genrotangle;
Leonardo Solis's avatar
Leonardo Solis committed
191
			}
192
			else	// If rotating around rotatable bond
Leonardo Solis's avatar
Leonardo Solis committed
193
			{
194
				uint rotbond_id = (rotation_list_element & RLIST_RBONDID_MASK) >> RLIST_RBONDID_SHIFT;
Leonardo Solis's avatar
Leonardo Solis committed
195

lvs's avatar
lvs committed
196
197
198
				rotation_unitvec[0] = kerconst_conform->rotbonds_unit_vectors_const[3*rotbond_id];
				rotation_unitvec[1] = kerconst_conform->rotbonds_unit_vectors_const[3*rotbond_id+1];
				rotation_unitvec[2] = kerconst_conform->rotbonds_unit_vectors_const[3*rotbond_id+2];
199
				
lvs's avatar
lvs committed
200
201
202
				rotation_movingvec[0] = kerconst_conform->rotbonds_moving_vectors_const[3*rotbond_id];
				rotation_movingvec[1] = kerconst_conform->rotbonds_moving_vectors_const[3*rotbond_id+1];
				rotation_movingvec[2] = kerconst_conform->rotbonds_moving_vectors_const[3*rotbond_id+2];
Leonardo Solis's avatar
Leonardo Solis committed
203

204
205
				rotation_angle = genotype[6+rotbond_id]*DEG_TO_RAD;

206
207
				// Performing additionally the first movement which 
				// is needed only if rotating around rotatable bond
Leonardo Solis's avatar
Leonardo Solis committed
208
209
210
211
212
				atom_to_rotate[0] -= rotation_movingvec[0];
				atom_to_rotate[1] -= rotation_movingvec[1];
				atom_to_rotate[2] -= rotation_movingvec[2];
			}

213
214
215
216
217
218
219
220
			// Transforming orientation and torsion angles into quaternions
			rotation_angle  = rotation_angle * 0.5f;
			float sin_angle = native_sin(rotation_angle);
			quatrot_left_q  = native_cos(rotation_angle);
			quatrot_left_x  = sin_angle*rotation_unitvec[0];
			quatrot_left_y  = sin_angle*rotation_unitvec[1];
			quatrot_left_z  = sin_angle*rotation_unitvec[2];

221
222
223
224
			// Performing rotation
			if ((rotation_list_element & RLIST_GENROT_MASK) != 0)	// If general rotation,
										// two rotations should be performed
										// (multiplying the quaternions)
Leonardo Solis's avatar
Leonardo Solis committed
225
			{
226
227
				// Calculating quatrot_left*ref_orientation_quats_const,
				// which means that reference orientation rotation is the first
Leonardo Solis's avatar
Leonardo Solis committed
228
229
230
231
232
				quatrot_temp_q = quatrot_left_q;
				quatrot_temp_x = quatrot_left_x;
				quatrot_temp_y = quatrot_left_y;
				quatrot_temp_z = quatrot_left_z;

lvs's avatar
lvs committed
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
				quatrot_left_q = quatrot_temp_q*kerconst_conform->ref_orientation_quats_const[4*(*run_id)]-
						 quatrot_temp_x*kerconst_conform->ref_orientation_quats_const[4*(*run_id)+1]-
						 quatrot_temp_y*kerconst_conform->ref_orientation_quats_const[4*(*run_id)+2]-
						 quatrot_temp_z*kerconst_conform->ref_orientation_quats_const[4*(*run_id)+3];
				quatrot_left_x = quatrot_temp_q*kerconst_conform->ref_orientation_quats_const[4*(*run_id)+1]+
						 kerconst_conform->ref_orientation_quats_const[4*(*run_id)]*quatrot_temp_x+
						 quatrot_temp_y*kerconst_conform->ref_orientation_quats_const[4*(*run_id)+3]-
						 kerconst_conform->ref_orientation_quats_const[4*(*run_id)+2]*quatrot_temp_z;
				quatrot_left_y = quatrot_temp_q*kerconst_conform->ref_orientation_quats_const[4*(*run_id)+2]+
						 kerconst_conform->ref_orientation_quats_const[4*(*run_id)]*quatrot_temp_y+
						 kerconst_conform->ref_orientation_quats_const[4*(*run_id)+1]*quatrot_temp_z-
						 quatrot_temp_x*kerconst_conform->ref_orientation_quats_const[4*(*run_id)+3];
				quatrot_left_z = quatrot_temp_q*kerconst_conform->ref_orientation_quats_const[4*(*run_id)+3]+
						 kerconst_conform->ref_orientation_quats_const[4*(*run_id)]*quatrot_temp_z+
						 quatrot_temp_x*kerconst_conform->ref_orientation_quats_const[4*(*run_id)+2]-
						 kerconst_conform->ref_orientation_quats_const[4*(*run_id)+1]*quatrot_temp_y;
Leonardo Solis's avatar
Leonardo Solis committed
249
250
251
			}

			quatrot_temp_q = 0 -
252
253
254
					 quatrot_left_x*atom_to_rotate [0] -
					 quatrot_left_y*atom_to_rotate [1] -
					 quatrot_left_z*atom_to_rotate [2];
Leonardo Solis's avatar
Leonardo Solis committed
255
			quatrot_temp_x = quatrot_left_q*atom_to_rotate [0] +
256
257
					 quatrot_left_y*atom_to_rotate [2] -
					 quatrot_left_z*atom_to_rotate [1];
Leonardo Solis's avatar
Leonardo Solis committed
258
			quatrot_temp_y = quatrot_left_q*atom_to_rotate [1] -
259
260
					 quatrot_left_x*atom_to_rotate [2] +
					 quatrot_left_z*atom_to_rotate [0];
Leonardo Solis's avatar
Leonardo Solis committed
261
			quatrot_temp_z = quatrot_left_q*atom_to_rotate [2] +
262
263
					 quatrot_left_x*atom_to_rotate [1] -
					 quatrot_left_y*atom_to_rotate [0];
Leonardo Solis's avatar
Leonardo Solis committed
264
265

			atom_to_rotate [0] = 0 -
266
267
268
269
					  quatrot_temp_q*quatrot_left_x +
					  quatrot_temp_x*quatrot_left_q -
					  quatrot_temp_y*quatrot_left_z +
					  quatrot_temp_z*quatrot_left_y;
Leonardo Solis's avatar
Leonardo Solis committed
270
			atom_to_rotate [1] = 0 -
271
272
273
274
					  quatrot_temp_q*quatrot_left_y +
					  quatrot_temp_x*quatrot_left_z +
					  quatrot_temp_y*quatrot_left_q -
					  quatrot_temp_z*quatrot_left_x;
Leonardo Solis's avatar
Leonardo Solis committed
275
			atom_to_rotate [2] = 0 -
276
277
278
279
					  quatrot_temp_q*quatrot_left_z -
					  quatrot_temp_x*quatrot_left_y +
					  quatrot_temp_y*quatrot_left_x +
					  quatrot_temp_z*quatrot_left_q;
Leonardo Solis's avatar
Leonardo Solis committed
280

281
			// Performing final movement and storing values
Leonardo Solis's avatar
Leonardo Solis committed
282
283
284
285
286
287
288
289
290
291
292
			calc_coords_x[atom_id] = atom_to_rotate [0] + rotation_movingvec[0];
			calc_coords_y[atom_id] = atom_to_rotate [1] + rotation_movingvec[1];
			calc_coords_z[atom_id] = atom_to_rotate [2] + rotation_movingvec[2];

		} // End if-statement not dummy rotation

		barrier(CLK_LOCAL_MEM_FENCE);

	} // End rotation_counter for-loop

	// ================================================
Leonardo Solis's avatar
Leonardo Solis committed
293
	// CALCULATING INTERMOLECULAR ENERGY
Leonardo Solis's avatar
Leonardo Solis committed
294
	// ================================================
295
296
297
	for (uint atom_id = get_local_id(0);
	          atom_id < dockpars_num_of_atoms;
	          atom_id+= NUM_OF_THREADS_PER_BLOCK)
Leonardo Solis's avatar
Leonardo Solis committed
298
	{
lvs's avatar
lvs committed
299
		uint atom_typeid = kerconst_interintra->atom_types_const[atom_id];
300
301
302
		float x = calc_coords_x[atom_id];
		float y = calc_coords_y[atom_id];
		float z = calc_coords_z[atom_id];
lvs's avatar
lvs committed
303
		float q = kerconst_interintra->atom_charges_const[atom_id];
Leonardo Solis's avatar
Leonardo Solis committed
304
305

		if ((x < 0) || (y < 0) || (z < 0) || (x >= dockpars_gridsize_x-1)
306
307
				                  || (y >= dockpars_gridsize_y-1)
						  || (z >= dockpars_gridsize_z-1)){
Leonardo Solis's avatar
Leonardo Solis committed
308
			partial_energies[get_local_id(0)] += 16777216.0f; //100000.0f;
309
	
lvs's avatar
lvs committed
310
			#if defined (DEBUG_ENERGY_KERNEL)
311
312
			partial_interE[get_local_id(0)] += 16777216.0f;
			#endif
Leonardo Solis's avatar
Leonardo Solis committed
313
314
315
		}
		else
		{
Leonardo Solis's avatar
Leonardo Solis committed
316
			// Getting coordinates
317
318
319
320
321
322
323
324
325
326
			int x_low  = (int)floor(x); 
			int y_low  = (int)floor(y); 
			int z_low  = (int)floor(z);
			int x_high = (int)ceil(x); 
			int y_high = (int)ceil(y); 
			int z_high = (int)ceil(z);
			float dx = x - x_low; 
			float dy = y - y_low; 
			float dz = z - z_low;

Leonardo Solis's avatar
Leonardo Solis committed
327
			// Calculating interpolation weights
328
			float weights[2][2][2];
Leonardo Solis's avatar
Leonardo Solis committed
329
330
331
332
333
334
335
336
337
			weights [0][0][0] = (1-dx)*(1-dy)*(1-dz);
			weights [1][0][0] = dx*(1-dy)*(1-dz);
			weights [0][1][0] = (1-dx)*dy*(1-dz);
			weights [1][1][0] = dx*dy*(1-dz);
			weights [0][0][1] = (1-dx)*(1-dy)*dz;
			weights [1][0][1] = dx*(1-dy)*dz;
			weights [0][1][1] = (1-dx)*dy*dz;
			weights [1][1][1] = dx*dy*dz;

Leonardo Solis's avatar
Leonardo Solis committed
338
			// Capturing affinity values
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
			uint ylow_times_g1  = y_low*g1;
			uint yhigh_times_g1 = y_high*g1;
		  	uint zlow_times_g2  = z_low*g2;
			uint zhigh_times_g2 = z_high*g2;

			// Grid offset
			uint offset_cube_000 = x_low  + ylow_times_g1  + zlow_times_g2;
			uint offset_cube_100 = x_high + ylow_times_g1  + zlow_times_g2;
			uint offset_cube_010 = x_low  + yhigh_times_g1 + zlow_times_g2;
			uint offset_cube_110 = x_high + yhigh_times_g1 + zlow_times_g2;
			uint offset_cube_001 = x_low  + ylow_times_g1  + zhigh_times_g2;
			uint offset_cube_101 = x_high + ylow_times_g1  + zhigh_times_g2;
			uint offset_cube_011 = x_low  + yhigh_times_g1 + zhigh_times_g2;
			uint offset_cube_111 = x_high + yhigh_times_g1 + zhigh_times_g2;

			uint mul_tmp = atom_typeid*g3;

			float cube[2][2][2];
			cube [0][0][0] = *(dockpars_fgrids + offset_cube_000 + mul_tmp);
			cube [1][0][0] = *(dockpars_fgrids + offset_cube_100 + mul_tmp);
			cube [0][1][0] = *(dockpars_fgrids + offset_cube_010 + mul_tmp);
		        cube [1][1][0] = *(dockpars_fgrids + offset_cube_110 + mul_tmp);
		        cube [0][0][1] = *(dockpars_fgrids + offset_cube_001 + mul_tmp);
			cube [1][0][1] = *(dockpars_fgrids + offset_cube_101 + mul_tmp);
                        cube [0][1][1] = *(dockpars_fgrids + offset_cube_011 + mul_tmp);
                        cube [1][1][1] = *(dockpars_fgrids + offset_cube_111 + mul_tmp);
Leonardo Solis's avatar
Leonardo Solis committed
365

Leonardo Solis's avatar
Leonardo Solis committed
366
			// Calculating affinity energy
Leonardo Solis's avatar
Leonardo Solis committed
367
368
			partial_energies[get_local_id(0)] += TRILININTERPOL(cube, weights);

lvs's avatar
lvs committed
369
			#if defined (DEBUG_ENERGY_KERNEL)
370
371
372
			partial_interE[get_local_id(0)] += TRILININTERPOL(cube, weights);
			#endif

Leonardo Solis's avatar
Leonardo Solis committed
373
			// Capturing electrostatic values
374
			atom_typeid = dockpars_num_of_atypes;
Leonardo Solis's avatar
Leonardo Solis committed
375

376
377
378
379
380
381
382
383
384
			mul_tmp = atom_typeid*g3;
			cube [0][0][0] = *(dockpars_fgrids + offset_cube_000 + mul_tmp);
			cube [1][0][0] = *(dockpars_fgrids + offset_cube_100 + mul_tmp);
      			cube [0][1][0] = *(dockpars_fgrids + offset_cube_010 + mul_tmp);
      			cube [1][1][0] = *(dockpars_fgrids + offset_cube_110 + mul_tmp);
		       	cube [0][0][1] = *(dockpars_fgrids + offset_cube_001 + mul_tmp);
		        cube [1][0][1] = *(dockpars_fgrids + offset_cube_101 + mul_tmp);
		        cube [0][1][1] = *(dockpars_fgrids + offset_cube_011 + mul_tmp);
		        cube [1][1][1] = *(dockpars_fgrids + offset_cube_111 + mul_tmp);
385
386

			// Calculating electrostatic energy
Leonardo Solis's avatar
Leonardo Solis committed
387
388
			partial_energies[get_local_id(0)] += q * TRILININTERPOL(cube, weights);

lvs's avatar
lvs committed
389
			#if defined (DEBUG_ENERGY_KERNEL)
390
391
392
			partial_interE[get_local_id(0)] += q * TRILININTERPOL(cube, weights);
			#endif

Leonardo Solis's avatar
Leonardo Solis committed
393
			// Capturing desolvation values
394
			atom_typeid = dockpars_num_of_atypes+1;
Leonardo Solis's avatar
Leonardo Solis committed
395

396
397
398
399
400
401
402
403
404
			mul_tmp = atom_typeid*g3;
			cube [0][0][0] = *(dockpars_fgrids + offset_cube_000 + mul_tmp);
			cube [1][0][0] = *(dockpars_fgrids + offset_cube_100 + mul_tmp);
      			cube [0][1][0] = *(dockpars_fgrids + offset_cube_010 + mul_tmp);
      			cube [1][1][0] = *(dockpars_fgrids + offset_cube_110 + mul_tmp);
      			cube [0][0][1] = *(dockpars_fgrids + offset_cube_001 + mul_tmp);
      			cube [1][0][1] = *(dockpars_fgrids + offset_cube_101 + mul_tmp);
      			cube [0][1][1] = *(dockpars_fgrids + offset_cube_011 + mul_tmp);
      			cube [1][1][1] = *(dockpars_fgrids + offset_cube_111 + mul_tmp);
Leonardo Solis's avatar
Leonardo Solis committed
405

Leonardo Solis's avatar
Leonardo Solis committed
406
			// Calculating desolvation energy
Leonardo Solis's avatar
Leonardo Solis committed
407
			partial_energies[get_local_id(0)] += fabs(q) * TRILININTERPOL(cube, weights);
408

lvs's avatar
lvs committed
409
			#if defined (DEBUG_ENERGY_KERNEL)
410
411
			partial_interE[get_local_id(0)] += fabs(q) * TRILININTERPOL(cube, weights);
			#endif
Leonardo Solis's avatar
Leonardo Solis committed
412
413
		}

414
	} // End atom_id for-loop (INTERMOLECULAR ENERGY)
Leonardo Solis's avatar
Leonardo Solis committed
415

416

lvs's avatar
lvs committed
417
	#if defined (DEBUG_ENERGY_KERNEL)
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
	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


Leonardo Solis's avatar
Leonardo Solis committed
438
439
	// In paper: intermolecular and internal energy calculation
	// are independent from each other, -> NO BARRIER NEEDED
440
  	// but require different operations,
Leonardo Solis's avatar
Leonardo Solis committed
441
442
443
	// thus, they can be executed only sequentially on the GPU.

	// ================================================
Leonardo Solis's avatar
Leonardo Solis committed
444
	// CALCULATING INTRAMOLECULAR ENERGY
Leonardo Solis's avatar
Leonardo Solis committed
445
	// ================================================
446
447
448
	for (uint contributor_counter = get_local_id(0);
	          contributor_counter < dockpars_num_of_intraE_contributors;
	          contributor_counter +=NUM_OF_THREADS_PER_BLOCK)
lvs's avatar
lvs committed
449
450
451
452
453
454
455

#if 0
if (get_local_id (0) == 0) {
	for (uint contributor_counter = 0;
	          contributor_counter < dockpars_num_of_intraE_contributors;
	          contributor_counter ++)
#endif
Leonardo Solis's avatar
Leonardo Solis committed
456
	{
lvs's avatar
lvs committed
457
458
459
460
461
462
463
464
465
466
#if 0
		// Only for testing smoothing
		float smoothed_intraE = 0.0f;

		float raw_intraE_vdw_hb = 0.0f;
		float raw_intraE_el     = 0.0f;
		float raw_intraE_sol    = 0.0f;
		float raw_intraE        = 0.0f;
#endif

Leonardo Solis's avatar
Leonardo Solis committed
467
		// Getting atom IDs
lvs's avatar
lvs committed
468
469
		uint atom1_id = kerconst_intracontrib->intraE_contributors_const[3*contributor_counter];
		uint atom2_id = kerconst_intracontrib->intraE_contributors_const[3*contributor_counter+1];
Leonardo Solis's avatar
Leonardo Solis committed
470

471
472
473
474
475
		// Calculating vector components of vector going
		// from first atom's to second atom's coordinates
		float subx = calc_coords_x[atom1_id] - calc_coords_x[atom2_id];
		float suby = calc_coords_y[atom1_id] - calc_coords_y[atom2_id];
		float subz = calc_coords_z[atom1_id] - calc_coords_z[atom2_id];
Leonardo Solis's avatar
Leonardo Solis committed
476

Leonardo Solis's avatar
Leonardo Solis committed
477
		// Calculating atomic_distance
478
		float atomic_distance = native_sqrt(subx*subx + suby*suby + subz*subz)*dockpars_grid_spacing;
Leonardo Solis's avatar
Leonardo Solis committed
479

lvs's avatar
lvs committed
480
		// Getting type IDs
lvs's avatar
lvs committed
481
482
		uint atom1_typeid = kerconst_interintra->atom_types_const[atom1_id];
		uint atom2_typeid = kerconst_interintra->atom_types_const[atom2_id];
Leonardo Solis's avatar
Leonardo Solis committed
483

lvs's avatar
lvs committed
484
485
		uint atom1_type_vdw_hb = kerconst_intra->atom1_types_reqm_const [atom1_typeid];
     	        uint atom2_type_vdw_hb = kerconst_intra->atom2_types_reqm_const [atom2_typeid];
lvs's avatar
lvs committed
486

lvs's avatar
lvs committed
487
488
489
490
491
492
		// 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;
lvs's avatar
lvs committed
493

lvs's avatar
lvs committed
494
		if (kerconst_intracontrib->intraE_contributors_const[3*contributor_counter+2] == 1)	//H-bond
lvs's avatar
lvs committed
495
		{
lvs's avatar
lvs committed
496
			opt_distance = kerconst_intra->reqm_hbond_const [atom1_type_vdw_hb] + kerconst_intra->reqm_hbond_const [atom2_type_vdw_hb];
lvs's avatar
lvs committed
497
498
499
		}
		else	//van der Waals
		{
lvs's avatar
lvs committed
500
			opt_distance = 0.5f*(kerconst_intra->reqm_const [atom1_type_vdw_hb] + kerconst_intra->reqm_const [atom2_type_vdw_hb]);
lvs's avatar
lvs committed
501
		}
lvs's avatar
lvs committed
502

lvs's avatar
lvs committed
503
504
505
		// Getting smoothed distance
		// smoothed_distance = function(atomic_distance, opt_distance)
		float smoothed_distance;
506
		float delta_distance = 0.5f*dockpars_smooth; 
lvs's avatar
lvs committed
507

lvs's avatar
lvs committed
508
509
510
511
512
513
514
515
516
		if (atomic_distance <= (opt_distance - delta_distance)) {
			smoothed_distance = atomic_distance + delta_distance;
		}
		else if (atomic_distance < (opt_distance + delta_distance)) {
			smoothed_distance = opt_distance;
		}
		else { // else if (atomic_distance >= (opt_distance + delta_distance))
			smoothed_distance = atomic_distance - delta_distance;
		}
lvs's avatar
lvs committed
517

lvs's avatar
lvs committed
518
		// Calculating energy contributions
519
		// Cuttoff1: internuclear-distance at 8A only for vdw and hbond.
lvs's avatar
lvs committed
520
521
		if (atomic_distance < 8.0f)
		{
Leonardo Solis's avatar
Leonardo Solis committed
522
			// Calculating van der Waals / hydrogen bond term
lvs's avatar
lvs committed
523
			partial_energies[get_local_id(0)] += native_divide(kerconst_intra->VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance/*atomic_distance*/,12));
524
525

			#if 0
lvs's avatar
lvs committed
526
527
			smoothed_intraE = native_divide(kerconst_intra->VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance,12));
			raw_intraE_vdw_hb      = native_divide(kerconst_intra->VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance  ,12)); 
528
529
			#endif

lvs's avatar
lvs committed
530
			#if defined (DEBUG_ENERGY_KERNEL)
lvs's avatar
lvs committed
531
			partial_intraE[get_local_id(0)] += native_divide(kerconst_intra->VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance/*atomic_distance*/,12));
532
533
			#endif

lvs's avatar
lvs committed
534
535
			if (kerconst_intracontrib->intraE_contributors_const[3*contributor_counter+2] == 1) {	//H-bond
				partial_energies[get_local_id(0)] -= native_divide(kerconst_intra->VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance/*atomic_distance*/,10));
536
537

				#if 0
lvs's avatar
lvs committed
538
539
				smoothed_intraE -= native_divide(kerconst_intra->VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance,10));
				raw_intraE_vdw_hb 	-= native_divide(kerconst_intra->VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance  ,10));
540
541
				#endif

lvs's avatar
lvs committed
542
				#if defined (DEBUG_ENERGY_KERNEL)
lvs's avatar
lvs committed
543
				partial_intraE[get_local_id(0)] -= native_divide(kerconst_intra->VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance/*atomic_distance*/,10));
544
				#endif
545
546
			}
			else {	//van der Waals
lvs's avatar
lvs committed
547
				partial_energies[get_local_id(0)] -= native_divide(kerconst_intra->VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance/*atomic_distance*/,6));
548
549

				#if 0
lvs's avatar
lvs committed
550
551
				smoothed_intraE -= native_divide(kerconst_intra->VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance,6));
				raw_intraE_vdw_hb      -= native_divide(kerconst_intra->VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance  ,6));
552
553
				#endif

lvs's avatar
lvs committed
554
				#if defined (DEBUG_ENERGY_KERNEL)
lvs's avatar
lvs committed
555
				partial_intraE[get_local_id(0)] -= native_divide(kerconst_intra->VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance/*atomic_distance*/,6));
556
				#endif
557
			}
558
		} // if cuttoff1 - internuclear-distance at 8A	
lvs's avatar
lvs committed
559

560
561
562
563
564
565
		// Calculating energy contributions
		// Cuttoff2: internuclear-distance at 20.48A only for el and sol.
		if (atomic_distance < 20.48f)
		{
			// Calculating electrostatic term
	       		partial_energies[get_local_id(0)] += native_divide (
lvs's avatar
lvs committed
566
		                                                     dockpars_coeff_elec * kerconst_interintra->atom_charges_const[atom1_id] * kerconst_interintra->atom_charges_const[atom2_id],
567
568
569
570
		                                                     atomic_distance * (DIEL_A + native_divide(DIEL_B,(1.0f + DIEL_K*native_exp(-DIEL_B_TIMES_H*atomic_distance))))
		                                                     );
			#if 0
			smoothed_intraE += native_divide (
lvs's avatar
lvs committed
571
						      dockpars_coeff_elec * kerconst_interintra->atom_charges_const[atom1_id] * kerconst_interintra->atom_charges_const[atom2_id],
572
573
574
575
						      atomic_distance * (DIEL_A + native_divide(DIEL_B,(1.0f + DIEL_K*native_exp(-DIEL_B_TIMES_H*atomic_distance))))
						     );

			raw_intraE_el 	= native_divide (
lvs's avatar
lvs committed
576
						      dockpars_coeff_elec * kerconst_interintra->atom_charges_const[atom1_id] * kerconst_interintra->atom_charges_const[atom2_id],
577
578
579
580
581
582
						      atomic_distance * (DIEL_A + native_divide(DIEL_B,(1.0f + DIEL_K*native_exp(-DIEL_B_TIMES_H*atomic_distance))))
						     );
			#endif

			#if defined (DEBUG_ENERGY_KERNEL)
			partial_intraE[get_local_id(0)] += native_divide (
583
		                                                     dockpars_coeff_elec * kerconst_interintra->atom_charges_const[atom1_id] * kerconst_interintra->atom_charges_const[atom2_id],
584
585
586
		                                                     atomic_distance * (DIEL_A + native_divide(DIEL_B,(1.0f + DIEL_K*native_exp(-DIEL_B_TIMES_H*atomic_distance))))
		                                                     );
			#endif
lvs's avatar
lvs committed
587

588
589
			// Calculating desolvation term
			// 1/25.92 = 0.038580246913580245
lvs's avatar
lvs committed
590
591
592
593
			partial_energies[get_local_id(0)] += ((kerconst_intra->dspars_S_const[atom1_typeid] +
								       dockpars_qasp*fabs(kerconst_interintra->atom_charges_const[atom1_id]))*kerconst_intra->dspars_V_const[atom2_typeid] +
							               (kerconst_intra->dspars_S_const[atom2_typeid] +
								       dockpars_qasp*fabs(kerconst_interintra->atom_charges_const[atom2_id]))*kerconst_intra->dspars_V_const[atom1_typeid]) *
594
595
596
							               dockpars_coeff_desolv*native_exp(-0.03858025f*native_powr(atomic_distance, 2));

			#if 0
lvs's avatar
lvs committed
597
598
599
600
			smoothed_intraE += ((kerconst_intra->dspars_S_const[atom1_typeid] +
					dockpars_qasp*fabs(kerconst_interintra->atom_charges_const[atom1_id]))*kerconst_intra->dspars_V_const[atom2_typeid] +
				       (kerconst_intra->dspars_S_const[atom2_typeid] +
				       dockpars_qasp*fabs(kerconst_interintra->atom_charges_const[atom2_id]))*kerconst_intra->dspars_V_const[atom1_typeid]) *
601
602
				       dockpars_coeff_desolv*native_exp(-0.03858025f*native_powr(atomic_distance, 2));

lvs's avatar
lvs committed
603
604
605
606
			raw_intraE_sol = ((kerconst_intra->dspars_S_const[atom1_typeid] +
					dockpars_qasp*fabs(kerconst_interintra->atom_charges_const[atom1_id]))*kerconst_intra->dspars_V_const[atom2_typeid] +
				       (kerconst_intra->dspars_S_const[atom2_typeid] +
				       dockpars_qasp*fabs(kerconst_interintra->atom_charges_const[atom2_id]))*kerconst_intra->dspars_V_const[atom1_typeid]) *
607
608
				       dockpars_coeff_desolv*native_exp(-0.03858025f*native_powr(atomic_distance, 2));
			#endif
609

610
			#if defined (DEBUG_ENERGY_KERNEL)
lvs's avatar
lvs committed
611
612
613
614
			partial_intraE[get_local_id(0)] += ((kerconst_intra->dspars_S_const[atom1_typeid] +
							       dockpars_qasp*fabs(kerconst_interintra->atom_charges_const[atom1_id]))*kerconst_intra->dspars_V_const[atom2_typeid] +
						               (kerconst_intra->dspars_S_const[atom2_typeid] +
							       dockpars_qasp*fabs(kerconst_interintra->atom_charges_const[atom2_id]))*kerconst_intra->dspars_V_const[atom1_typeid]) *
615
616
617
						               dockpars_coeff_desolv*native_exp(-0.03858025f*native_powr(atomic_distance, 2));
			#endif
		} // if cuttoff2 - internuclear-distance at 20.48A
lvs's avatar
lvs committed
618
619
620
621
622
623
624
625
626
627



#if 0
			raw_intraE = raw_intraE_vdw_hb + raw_intraE_el + raw_intraE_sol;

			if (debug == true) {
///*
			//if (get_local_id (0) == 0) {

lvs's avatar
lvs committed
628
				if (kerconst_intracontrib->intraE_contributors_const[3*contributor_counter+2] == 1)	//H-bond
lvs's avatar
lvs committed
629
630
				{
					// diogos table
lvs's avatar
lvs committed
631
//					printf("%3u %-5s %3u %3u %3.5f %3.5f %3.5f %3.5f %3.5f %3.5f %3.5f\n", contributor_counter, "hbond", atom1_id, atom2_id, kerconst_intra->reqm_hbond_const [atom1_type_vdw_hb], kerconst_intra->reqm_hbond_const [atom2_type_vdw_hb], opt_distance, atomic_distance, smoothed_distance, smoothed_intraE, raw_intraE);
lvs's avatar
lvs committed
632

lvs's avatar
lvs committed
633
					printf("%3u %-5s %3u %3u %3.5f %3.5f %3.5f %3.5f %3.5f %3.5f %3.5f %3.5f %3.5f %3.5f\n", contributor_counter, "hbond", atom1_id, atom2_id, kerconst_intra->reqm_hbond_const [atom1_type_vdw_hb], kerconst_intra->reqm_hbond_const [atom2_type_vdw_hb], opt_distance, atomic_distance, smoothed_distance, smoothed_intraE, raw_intraE, raw_intraE_vdw_hb, raw_intraE_el, raw_intraE_sol);
lvs's avatar
lvs committed
634
635
636
637
638
639
640
641




				}
				else	//van der Waals
				{
					// diogos table
lvs's avatar
lvs committed
642
//					printf("%3u %-5s %3u %3u %3.5f %3.5f %3.5f %3.5f %3.5f %3.5f %3.5f\n", contributor_counter, "vdw", atom1_id, atom2_id, kerconst_intra->reqm_const [atom1_type_vdw_hb], kerconst_intra->reqm_const [atom2_type_vdw_hb], opt_distance, atomic_distance, smoothed_distance, smoothed_intraE, raw_intraE);
lvs's avatar
lvs committed
643

lvs's avatar
lvs committed
644
					printf("%3u %-5s %3u %3u %3.5f %3.5f %3.5f %3.5f %3.5f %3.5f %3.5f %3.5f %3.5f %3.5f\n", contributor_counter, "vdw", atom1_id, atom2_id, kerconst_intra->reqm_const [atom1_type_vdw_hb], kerconst_intra->reqm_const [atom2_type_vdw_hb], opt_distance, atomic_distance, smoothed_distance, smoothed_intraE, raw_intraE, raw_intraE_vdw_hb, raw_intraE_el, raw_intraE_sol);
lvs's avatar
lvs committed
645
646
647
648
649
650
651
652


				}; 
			//}
//*/
			}
#endif

653
654
	} // End contributor_counter for-loop (INTRAMOLECULAR ENERGY)

lvs's avatar
lvs committed
655
656
657
658
#if 0
} // if (get_local_id (0) == 0) {
#endif

Leonardo Solis's avatar
Leonardo Solis committed
659
660
661
662
663
664
	barrier(CLK_LOCAL_MEM_FENCE);

	if (get_local_id(0) == 0)
	{
		*energy = partial_energies[0];

665
666
667
		for (uint contributor_counter=1;
		          contributor_counter<NUM_OF_THREADS_PER_BLOCK;
		          contributor_counter++)
Leonardo Solis's avatar
Leonardo Solis committed
668
669
670
671
		{
			*energy += partial_energies[contributor_counter];
		}
	}
672

673
674
	barrier(CLK_LOCAL_MEM_FENCE);

lvs's avatar
lvs committed
675
	#if defined (DEBUG_ENERGY_KERNEL)
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
	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

Leonardo Solis's avatar
Leonardo Solis committed
693
694
}

695
696
// No needed to be included as all kernel sources are stringified
#if 0
Leonardo Solis's avatar
Leonardo Solis committed
697
698
699
700
701
#include "kernel1.cl"
#include "kernel2.cl"
#include "auxiliary_genetic.cl"
#include "kernel4.cl"
#include "kernel3.cl"
702
#include "calcgradient.cl"
lvs's avatar
lvs committed
703
#include "kernel_sd.cl"
lvs's avatar
lvs committed
704
#include "kernel_fire.cl"
lvs's avatar
lvs committed
705
#include "kernel_ad.cl"
706
707
708
#include "calcEnerGrad.cl"
#endif