calcenergy.cl 32.3 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
24
/*

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.

*/


25
26
27
#define DEBUG_ENERGY


28
#include "calcenergy_basic.h"
Leonardo Solis's avatar
Leonardo Solis committed
29
30
31

// All related pragmas are in defines.h (accesible by host and device code)

32
33
34
35
36
37
38
39
40
41
42
43
44
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,
		 __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,
Leonardo Solis's avatar
Leonardo Solis committed
45

46
47
48
49
                    // 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.
50
51
52
53
54
55
56
57
58
		    	__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,

59
60
61
62
63
			#if defined (DEBUG_ENERGY)
			__local float* partial_interE,
			__local float* partial_intraE,
			#endif

64
65
66
67
68
69
70
71
72
73
74
75
76
	             __constant float* atom_charges_const,
                     __constant char*  atom_types_const,
                     __constant char*  intraE_contributors_const,
                     __constant float* VWpars_AC_const,
                     __constant float* VWpars_BD_const,
                     __constant float* dspars_S_const,
                     __constant float* dspars_V_const,
                     __constant int*   rotlist_const,
                     __constant float* ref_coords_x_const,
                     __constant float* ref_coords_y_const,
                     __constant float* ref_coords_z_const,
                     __constant float* rotbonds_moving_vectors_const,
                     __constant float* rotbonds_unit_vectors_const,
77
                     __constant float* ref_orientation_quats_const
Leonardo Solis's avatar
Leonardo Solis committed
78
79
80
81
82
83
84
85
86
)

//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;

87
88
89
90
91
92
	#if defined (DEBUG_ENERGY)
	partial_interE[get_local_id(0)] = 0.0f;
	partial_intraE[get_local_id(0)] = 0.0f;
	#endif


93
#if defined (IMPROVE_GRID)
Leonardo Solis's avatar
Leonardo Solis committed
94
95
	// INTERMOLECULAR for-loop (intermediate results)
	// It stores a product of two chars
96
	//uint mul_tmp;
Leonardo Solis's avatar
Leonardo Solis committed
97

98
99
100
	uchar g1 = dockpars_gridsize_x;
	uint  g2 = dockpars_gridsize_x * dockpars_gridsize_y;
  	uint  g3 = dockpars_gridsize_x * dockpars_gridsize_y * dockpars_gridsize_z;
Leonardo Solis's avatar
Leonardo Solis committed
101
#else
102

103
#endif
104

Leonardo Solis's avatar
Leonardo Solis committed
105
	// ================================================
Leonardo Solis's avatar
Leonardo Solis committed
106
	// CALCULATING ATOMIC POSITIONS AFTER ROTATIONS
Leonardo Solis's avatar
Leonardo Solis committed
107
	// ================================================
108
109
110
	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
111
	{
112
		int rotation_list_element = rotlist_const[rotation_counter];
Leonardo Solis's avatar
Leonardo Solis committed
113

114
		if ((rotation_list_element & RLIST_DUMMY_MASK) == 0)	// If not dummy rotation
Leonardo Solis's avatar
Leonardo Solis committed
115
		{
116
117
118
119
			uint atom_id = rotation_list_element & RLIST_ATOMID_MASK;

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

121
			if ((rotation_list_element & RLIST_FIRSTROT_MASK) != 0)	// If first rotation of this atom
Leonardo Solis's avatar
Leonardo Solis committed
122
123
124
125
126
127
128
129
130
131
132
133
			{
				atom_to_rotate[0] = ref_coords_x_const[atom_id];
				atom_to_rotate[1] = ref_coords_y_const[atom_id];
				atom_to_rotate[2] = ref_coords_z_const[atom_id];
			}
			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];
			}

134
135
136
137
138
139
140
			// Capturing rotation vectors and angle
			float rotation_movingvec[3];

			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
141
			{
142
143
144
				// -------------------------------------------------------------------
				// Replacing rotation genes: from spherical space to Shoemake space
				// gene [0:2]: translation -> kept as original x, y, z
Leonardo Solis's avatar
Leonardo Solis committed
145
				// gene [3:5]: rotation    -> transformed into Shoemake (u1, u2, u3)
146
				// gene [6:N]: torsions	   -> kept as original angles	(all in sexagesimal)
147

148
				// Shoemake ranges:
Leonardo Solis's avatar
Leonardo Solis committed
149
				// u1, u2, u3: [0, 1]
150

151
152
153
				// Random generator in the host is changed:
				// LCG (original, myrand()) -> CPP std (rand())
				// -------------------------------------------------------------------
154

155
156
157
				// Transforming Shoemake (u1, u2, u3) genes into quaternions
				// FIXME: add precision choices with preprocessor directives: 
				// NATIVE_PRECISION, HALF_PRECISION, Full precision
158

159
				// Rotational genes in the Shoemake space expressed in radians
160
				float u1 = genotype[3];
Leonardo Solis's avatar
Leonardo Solis committed
161
162
				float u2 = genotype[4];
				float u3 = genotype[5];
163
164
165
166
167
168

				// u1, u2, u3 should be within their valid range of [0,1]
				quatrot_left_q = native_sqrt(1 - u1) * native_sin(PI_TIMES_2*u2); 
				quatrot_left_x = native_sqrt(1 - u1) * native_cos(PI_TIMES_2*u2);
				quatrot_left_y = native_sqrt(u1)     * native_sin(PI_TIMES_2*u3);
				quatrot_left_z = native_sqrt(u1)     * native_cos(PI_TIMES_2*u3);
169

Leonardo Solis's avatar
Leonardo Solis committed
170
171
172
173
				rotation_movingvec[0] = genotype[0];
				rotation_movingvec[1] = genotype[1];
				rotation_movingvec[2] = genotype[2];
			}
174
			else	// If rotating around rotatable bond
Leonardo Solis's avatar
Leonardo Solis committed
175
			{
176
				uint rotbond_id = (rotation_list_element & RLIST_RBONDID_MASK) >> RLIST_RBONDID_SHIFT;
Leonardo Solis's avatar
Leonardo Solis committed
177

178
				float rotation_unitvec[3];
Leonardo Solis's avatar
Leonardo Solis committed
179
180
181
				rotation_unitvec[0] = rotbonds_unit_vectors_const[3*rotbond_id];
				rotation_unitvec[1] = rotbonds_unit_vectors_const[3*rotbond_id+1];
				rotation_unitvec[2] = rotbonds_unit_vectors_const[3*rotbond_id+2];
182
				float rotation_angle = genotype[6+rotbond_id]*DEG_TO_RAD;
Leonardo Solis's avatar
Leonardo Solis committed
183
184
185
186
187

				rotation_movingvec[0] = rotbonds_moving_vectors_const[3*rotbond_id];
				rotation_movingvec[1] = rotbonds_moving_vectors_const[3*rotbond_id+1];
				rotation_movingvec[2] = rotbonds_moving_vectors_const[3*rotbond_id+2];

188
189
				// Performing additionally the first movement which 
				// is needed only if rotating around rotatable bond
Leonardo Solis's avatar
Leonardo Solis committed
190
191
192
				atom_to_rotate[0] -= rotation_movingvec[0];
				atom_to_rotate[1] -= rotation_movingvec[1];
				atom_to_rotate[2] -= rotation_movingvec[2];
193
194
195
196

				// Transforming torsion angles into quaternions
				// FIXME: add precision choices with preprocessor directives: 
				// NATIVE_PRECISION, HALF_PRECISION, Full precision
197
198
199
200
201
202
				rotation_angle  = native_divide(rotation_angle, 2.0f);
				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];
Leonardo Solis's avatar
Leonardo Solis committed
203
204
			}

205
206
207
208
			// 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
209
			{
210
211
				// Calculating quatrot_left*ref_orientation_quats_const,
				// which means that reference orientation rotation is the first
Leonardo Solis's avatar
Leonardo Solis committed
212
213
214
215
216
217
				quatrot_temp_q = quatrot_left_q;
				quatrot_temp_x = quatrot_left_x;
				quatrot_temp_y = quatrot_left_y;
				quatrot_temp_z = quatrot_left_z;

				quatrot_left_q = quatrot_temp_q*ref_orientation_quats_const[4*(*run_id)]-
218
219
220
						 quatrot_temp_x*ref_orientation_quats_const[4*(*run_id)+1]-
						 quatrot_temp_y*ref_orientation_quats_const[4*(*run_id)+2]-
						 quatrot_temp_z*ref_orientation_quats_const[4*(*run_id)+3];
Leonardo Solis's avatar
Leonardo Solis committed
221
				quatrot_left_x = quatrot_temp_q*ref_orientation_quats_const[4*(*run_id)+1]+
222
223
224
						 ref_orientation_quats_const[4*(*run_id)]*quatrot_temp_x+
						 quatrot_temp_y*ref_orientation_quats_const[4*(*run_id)+3]-
						 ref_orientation_quats_const[4*(*run_id)+2]*quatrot_temp_z;
Leonardo Solis's avatar
Leonardo Solis committed
225
				quatrot_left_y = quatrot_temp_q*ref_orientation_quats_const[4*(*run_id)+2]+
226
227
228
						 ref_orientation_quats_const[4*(*run_id)]*quatrot_temp_y+
						 ref_orientation_quats_const[4*(*run_id)+1]*quatrot_temp_z-
						 quatrot_temp_x*ref_orientation_quats_const[4*(*run_id)+3];
Leonardo Solis's avatar
Leonardo Solis committed
229
				quatrot_left_z = quatrot_temp_q*ref_orientation_quats_const[4*(*run_id)+3]+
230
231
232
						 ref_orientation_quats_const[4*(*run_id)]*quatrot_temp_z+
						 quatrot_temp_x*ref_orientation_quats_const[4*(*run_id)+2]-
						 ref_orientation_quats_const[4*(*run_id)+1]*quatrot_temp_y;
Leonardo Solis's avatar
Leonardo Solis committed
233
234
235
			}

			quatrot_temp_q = 0 -
236
237
238
					 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
239
			quatrot_temp_x = quatrot_left_q*atom_to_rotate [0] +
240
241
					 quatrot_left_y*atom_to_rotate [2] -
					 quatrot_left_z*atom_to_rotate [1];
Leonardo Solis's avatar
Leonardo Solis committed
242
			quatrot_temp_y = quatrot_left_q*atom_to_rotate [1] -
243
244
					 quatrot_left_x*atom_to_rotate [2] +
					 quatrot_left_z*atom_to_rotate [0];
Leonardo Solis's avatar
Leonardo Solis committed
245
			quatrot_temp_z = quatrot_left_q*atom_to_rotate [2] +
246
247
					 quatrot_left_x*atom_to_rotate [1] -
					 quatrot_left_y*atom_to_rotate [0];
Leonardo Solis's avatar
Leonardo Solis committed
248
249

			atom_to_rotate [0] = 0 -
250
251
252
253
					  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
254
			atom_to_rotate [1] = 0 -
255
256
257
258
					  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
259
			atom_to_rotate [2] = 0 -
260
261
262
263
					  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
264

265
			// Performing final movement and storing values
Leonardo Solis's avatar
Leonardo Solis committed
266
267
268
269
270
271
272
273
274
275
276
			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
277
	// CALCULATING INTERMOLECULAR ENERGY
Leonardo Solis's avatar
Leonardo Solis committed
278
	// ================================================
279
280
281
	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
282
	{
283
284
285
286
287
		uint atom_typeid = atom_types_const[atom_id];
		float x = calc_coords_x[atom_id];
		float y = calc_coords_y[atom_id];
		float z = calc_coords_z[atom_id];
		float q = atom_charges_const[atom_id];
Leonardo Solis's avatar
Leonardo Solis committed
288
289

		if ((x < 0) || (y < 0) || (z < 0) || (x >= dockpars_gridsize_x-1)
290
291
				                  || (y >= dockpars_gridsize_y-1)
						  || (z >= dockpars_gridsize_z-1)){
Leonardo Solis's avatar
Leonardo Solis committed
292
			partial_energies[get_local_id(0)] += 16777216.0f; //100000.0f;
293
294
295
296
	
			#if defined (DEBUG_ENERGY)
			partial_interE[get_local_id(0)] += 16777216.0f;
			#endif
Leonardo Solis's avatar
Leonardo Solis committed
297
298
299
		}
		else
		{
Leonardo Solis's avatar
Leonardo Solis committed
300
			// Getting coordinates
301
302
303
304
305
306
307
308
309
310
			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
311
			// Calculating interpolation weights
312
			float weights[2][2][2];
Leonardo Solis's avatar
Leonardo Solis committed
313
314
315
316
317
318
319
320
321
			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
322
			// Capturing affinity values
Leonardo Solis's avatar
Leonardo Solis committed
323
#if defined (IMPROVE_GRID)
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
			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
350
#else
351
352
353
354
355
356
357
			// -------------------------------------------------------------------
			// FIXME: this block within the "#else" preprocessor directive 
			// provides NO gradient corresponding to "atype" intermolecular energy
			// -------------------------------------------------------------------	

			cube [0][0][0] = GETGRIDVALUE(dockpars_fgrids, 
						      dockpars_gridsize_x,
358
						      dockparsdockpars_num_of_atoms;_gridsize_y, 
359
						      dockpars_gridsize_z,
360
						      atom_typeid, z_low, y_low, x_low);
361
362
363
364
			cube [1][0][0] = GETGRIDVALUE(dockpars_fgrids, 
						      dockpars_gridsize_x,
						      dockpars_gridsize_y, 
						      dockpars_gridsize_z,
365
						      atom_typeid, z_low, y_low, x_high);
366
367
368
369
			cube [0][1][0] = GETGRIDVALUE(dockpars_fgrids, 
						      dockpars_gridsize_x,
						      dockpars_gridsize_y, 
						      dockpars_gridsize_z,
370
						      atom_typeid, z_low, y_high, x_low);
371
			cube [1][1][0] = GETGRIDVALUE(dockpars_fgrids, 
372
373
						      dockpars_gridsize_x,
						      dockpars_gridsize_y, 
374
						      dockpars_gridsize_z,
375
						      atom_typeid, z_low, y_high, x_high);
376
377
378
379
			cube [0][0][1] = GETGRIDVALUE(dockpars_fgrids, 
						      dockpars_gridsize_x,
						      dockpars_gridsize_y, 
						      dockpars_gridsize_z,
380
						      atom_typeid, z_high, y_low, x_low);
381
382
383
384
			cube [1][0][1] = GETGRIDVALUE(dockpars_fgrids, 
						      dockpars_gridsize_x,
						      dockpars_gridsize_y, 
						      dockpars_gridsize_z,
385
						      atom_typeid, z_high, y_low, x_high);
386
387
388
389
			cube [0][1][1] = GETGRIDVALUE(dockpars_fgrids, 
						      dockpars_gridsize_x,
						      dockpars_gridsize_y, 
						      dockpars_gridsize_z,
390
						      atom_typeid, z_high, y_high, x_low);
391
392
393
394
			cube [1][1][1] = GETGRIDVALUE(dockpars_fgrids, 
						      dockpars_gridsize_x,
						      dockpars_gridsize_y, 
						      dockpars_gridsize_z,
395
						      atom_typeid, z_high, y_high, x_high);
Leonardo Solis's avatar
Leonardo Solis committed
396
397
#endif

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

401
402
403
404
			#if defined (DEBUG_ENERGY)
			partial_interE[get_local_id(0)] += TRILININTERPOL(cube, weights);
			#endif

Leonardo Solis's avatar
Leonardo Solis committed
405
			// Capturing electrostatic values
406
			atom_typeid = dockpars_num_of_atypes;
Leonardo Solis's avatar
Leonardo Solis committed
407
408

#if defined (IMPROVE_GRID)
409
410
411
412
413
414
415
416
417
			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
418
#else
419
420
			// -------------------------------------------------------------------
			// FIXME: this block within the "#else" preprocessor directive 
421
			// provides NO gradient corresponding to "elec" intermolecular energybas
422
423
424
425
426
427
			// -------------------------------------------------------------------

			cube [0][0][0] = GETGRIDVALUE(dockpars_fgrids, 
						      dockpars_gridsize_x,
						      dockpars_gridsize_y, 
						      dockpars_gridsize_z,
428
					 	      atom_typeid, z_low, y_low, x_low);
429
430
431
432
			cube [1][0][0] = GETGRIDVALUE(dockpars_fgrids, 
						      dockpars_gridsize_x,
						      dockpars_gridsize_y,
						      dockpars_gridsize_z,
433
						      atom_typeid, z_low, y_low, x_high);
434
435
436
437
			cube [0][1][0] = GETGRIDVALUE(dockpars_fgrids,
						      dockpars_gridsize_x,
						      dockpars_gridsize_y,
                                                      dockpars_gridsize_z,
438
						      atom_typeid, z_low, y_high, x_low);
439
440
441
			cube [1][1][0] = GETGRIDVALUE(dockpars_fgrids,
						      dockpars_gridsize_x,
						      dockpars_gridsize_y,
442
						      dockpars_gridsize_z,bas
443
						      atom_typeid, z_low, y_high, x_high);
444
445
446
447
			cube [0][0][1] = GETGRIDVALUE(dockpars_fgrids,
						      dockpars_gridsize_x,
						      dockpars_gridsize_y,
						      dockpars_gridsize_z,
448
						      atom_typeid, z_high, y_low, x_low);
449
450
451
452
			cube [1][0][1] = GETGRIDVALUE(dockpars_fgrids,
						      dockpars_gridsize_x,
						      dockpars_gridsize_y,
						      dockpars_gridsize_z,
453
						      atom_typeid, z_high, y_low, x_high);
454
455
456
457
			cube [0][1][1] = GETGRIDVALUE(dockpars_fgrids,
						      dockpars_gridsize_x,
						      dockpars_gridsize_y,
						      dockpars_gridsize_z,
458
						      atom_typeid, z_high, y_high, x_low);
459
460
461
462
			cube [1][1][1] = GETGRIDVALUE(dockpars_fgrids,
						      dockpars_gridsize_x,
						      dockpars_gridsize_y,
						      dockpars_gridsize_z,
463
						      atom_typeid, z_high, y_high, x_high);
Leonardo Solis's avatar
Leonardo Solis committed
464
465
#endif

Leonardo Solis's avatar
Leonardo Solis committed
466
			// Calculating electrosatic energy
Leonardo Solis's avatar
Leonardo Solis committed
467
468
			partial_energies[get_local_id(0)] += q * TRILININTERPOL(cube, weights);

469
470
471
472
			#if defined (DEBUG_ENERGY)
			partial_interE[get_local_id(0)] += q * TRILININTERPOL(cube, weights);
			#endif

Leonardo Solis's avatar
Leonardo Solis committed
473
			// Capturing desolvation values
474
			atom_typeid = dockpars_num_of_atypes+1;
Leonardo Solis's avatar
Leonardo Solis committed
475
476

#if defined (IMPROVE_GRID)
477
478
479
480
481
482
483
484
485
			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
486
#else
487
488
489
490
491
492
493
494
495
			// -------------------------------------------------------------------
			// FIXME: this block within the "#else" preprocessor directive 
			// provides NO gradient corresponding to "dsol" intermolecular energy
			// -------------------------------------------------------------------

			cube [0][0][0] = GETGRIDVALUE(dockpars_fgrids,
						      dockpars_gridsize_x,
						      dockpars_gridsize_y,
						      dockpars_gridsize_z,
496
						      atom_typeid, z_low, y_low, x_low);
497
			cube [1][0][0] = GETGRIDVALUE(dockpars_fgridsu2,
498
499
500
						      dockpars_gridsize_x,
						      dockpars_gridsize_y,
						      dockpars_gridsize_z,
501
						      atom_typeid, z_low, y_low, x_high);
502
			cube [0][1][0] = GETGRIDVALUE(dockpars_fgrids,
503
						      dockpars_gridsize_x,rotbonds_unit_vectors_const
504
505
						      dockpars_gridsize_y,
						      dockpars_gridsize_z,
506
						      atom_typeid, z_low, y_high, x_low);
507
508
509
510
			cube [1][1][0] = GETGRIDVALUE(dockpars_fgrids,
						      dockpars_gridsize_x,
						      dockpars_gridsize_y,
						      dockpars_gridsize_z,
511
						      atom_typeid, z_low, y_high, x_high);
512
513
514
515
			cube [0][0][1] = GETGRIDVALUE(dockpars_fgrids, 
						      dockpars_gridsize_x,
						      dockpars_gridsize_y,
						      dockpars_gridsize_z,
516
						      atom_typeid, z_high, y_low, x_low);
517
518
519
520
			cube [1][0][1] = GETGRIDVALUE(dockpars_fgrids,
						      dockpars_gridsize_x,
						      dockpars_gridsize_y,
						      dockpars_gridsize_z,
521
						      atom_typeid, z_high, y_low, x_high);
522
523
524
525
			cube [0][1][1] = GETGRIDVALUE(dockpars_fgrids,
						      dockpars_gridsize_x,
						      dockpars_gridsize_y,
						      dockpars_gridsize_z,
526
						      atom_typeid, z_high, y_high, x_low);
527
528
529
530
			cube [1][1][1] = GETGRIDVALUE(dockpars_fgrids,
						      dockpars_gridsize_x,
						      dockpars_gridsize_y,
						      dockpars_gridsize_z,
531
						      atom_typeid, z_high, y_high, x_high);
Leonardo Solis's avatar
Leonardo Solis committed
532
533
#endif

Leonardo Solis's avatar
Leonardo Solis committed
534
			// Calculating desolvation energy
Leonardo Solis's avatar
Leonardo Solis committed
535
			partial_energies[get_local_id(0)] += fabs(q) * TRILININTERPOL(cube, weights);
536
537
538
539

			#if defined (DEBUG_ENERGY)
			partial_interE[get_local_id(0)] += fabs(q) * TRILININTERPOL(cube, weights);
			#endif
Leonardo Solis's avatar
Leonardo Solis committed
540
541
		}

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

544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565

	#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


Leonardo Solis's avatar
Leonardo Solis committed
566
567
	// In paper: intermolecular and internal energy calculation
	// are independent from each other, -> NO BARRIER NEEDED
568
  	// but require different operations,
Leonardo Solis's avatar
Leonardo Solis committed
569
570
571
	// thus, they can be executed only sequentially on the GPU.

	// ================================================
Leonardo Solis's avatar
Leonardo Solis committed
572
	// CALCULATING INTRAMOLECULAR ENERGY
Leonardo Solis's avatar
Leonardo Solis committed
573
	// ================================================
574
575
576
	for (uint contributor_counter = get_local_id(0);
	          contributor_counter < dockpars_num_of_intraE_contributors;
	          contributor_counter +=NUM_OF_THREADS_PER_BLOCK)
Leonardo Solis's avatar
Leonardo Solis committed
577
	{
Leonardo Solis's avatar
Leonardo Solis committed
578
		// Getting atom IDs
579
580
		uint atom1_id = intraE_contributors_const[3*contributor_counter];
		uint atom2_id = intraE_contributors_const[3*contributor_counter+1];
Leonardo Solis's avatar
Leonardo Solis committed
581

Leonardo Solis's avatar
Leonardo Solis committed
582
		// Calculating address of first atom's coordinates
583
584
585
		float subx = calc_coords_x[atom1_id];
		float suby = calc_coords_y[atom1_id];
		float subz = calc_coords_z[atom1_id];
Leonardo Solis's avatar
Leonardo Solis committed
586

Leonardo Solis's avatar
Leonardo Solis committed
587
		// Calculating address of second atom's coordinates
Leonardo Solis's avatar
Leonardo Solis committed
588
589
590
591
		subx -= calc_coords_x[atom2_id];
		suby -= calc_coords_y[atom2_id];
		subz -= calc_coords_z[atom2_id];

Leonardo Solis's avatar
Leonardo Solis committed
592
		// Calculating atomic_distance
Leonardo Solis's avatar
Leonardo Solis committed
593
#if defined (NATIVE_PRECISION)
594
		float atomic_distance = native_sqrt(subx*subx + suby*suby + subz*subz)*dockpars_grid_spacing;
Leonardo Solis's avatar
Leonardo Solis committed
595
#elif defined (HALF_PRECISION)
596
		float atomic_distance = half_sqrt(subx*subx + suby*suby + subz*subz)*dockpars_grid_spacing;
Leonardo Solis's avatar
Leonardo Solis committed
597
#else	// Full precision
598
		float atomic_distance = sqrt(subx*subx + suby*suby + subz*subz)*dockpars_grid_spacing;
Leonardo Solis's avatar
Leonardo Solis committed
599
600
#endif

601
602
		if (atomic_distance < 1.0f)
			atomic_distance = 1.0f;
Leonardo Solis's avatar
Leonardo Solis committed
603

Leonardo Solis's avatar
Leonardo Solis committed
604
		// Calculating energy contributions
605
		if ((atomic_distance < 8.0f) && (atomic_distance < 20.48f))
Leonardo Solis's avatar
Leonardo Solis committed
606
		{
607
608
609
			// Getting type IDs
			uint atom1_typeid = atom_types_const[atom1_id];
			uint atom2_typeid = atom_types_const[atom2_id];
Leonardo Solis's avatar
Leonardo Solis committed
610

Leonardo Solis's avatar
Leonardo Solis committed
611
			// Calculating van der Waals / hydrogen bond term
Leonardo Solis's avatar
Leonardo Solis committed
612
#if defined (NATIVE_PRECISION)
613
			partial_energies[get_local_id(0)] += native_divide(VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance,12));
614
615
616
617
618
619


			#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

Leonardo Solis's avatar
Leonardo Solis committed
620
#elif defined (HALF_PRECISION)
621
			for (uint rotation_counter = get_local_id(0);
622
623
	          rotation_counter < dockpars_rotbondlist_length;
	          rotation_counter+=NUM_OF_THREADS_PER_BLOCK)
624
	{partial_energies[get_local_id(0)] += half_divide(VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],half_powr(atomic_distance,12));
Leonardo Solis's avatar
Leonardo Solis committed
625
#else	// Full precision
626
			partial_energies[get_local_id(0)] += VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid]/powr(atomic_distance,12);
Leonardo Solis's avatar
Leonardo Solis committed
627
628
#endif

629
			if (intraE_contributors_const[3*contributor_counter+2] == 1) {	//H-bond
Leonardo Solis's avatar
Leonardo Solis committed
630
#if defined (NATIVE_PRECISION)
631
				partial_energies[get_local_id(0)] -= native_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance,10));
632
633
634
635
636
			#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


Leonardo Solis's avatar
Leonardo Solis committed
637
#elif defined (HALF_PRECISION)
638
				partial_energies[get_local_id(0)] -= half_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],half_powr(atomic_distance,10));
Leonardo Solis's avatar
Leonardo Solis committed
639
#else	// Full precision
640
				partial_energies[get_local_id(0)] -= VWpars_BD_const[atom1_typeid*dockpars_num_of_atypes+atom2_typeid]/powr(atomic_distance,10);
Leonardo Solis's avatar
Leonardo Solis committed
641
#endif
642
643
			}
			else {	//van der Waals
Leonardo Solis's avatar
Leonardo Solis committed
644
#if defined (NATIVE_PRECISION)
645
				partial_energies[get_local_id(0)] -= native_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance,6));
646
647
648
649
650

			#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

Leonardo Solis's avatar
Leonardo Solis committed
651
#elif defined (HALF_PRECISION)
652
				partial_energies[get_local_id(0)] -= half_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],half_powr(atomic_distance,6));
Leonardo Solis's avatar
Leonardo Solis committed
653
#else	// Full precision
654
				partial_energies[get_local_id(0)] -= VWpars_BD_const[atom1_typeid*dockpars_num_of_atypes+atom2_typeid]/powr(atomic_distance,6);
Leonardo Solis's avatar
Leonardo Solis committed
655
#endif
656
			}
Leonardo Solis's avatar
Leonardo Solis committed
657
			// Calculating electrostatic term
Leonardo Solis's avatar
Leonardo Solis committed
658
659

/*
Leonardo Solis's avatar
Leonardo Solis committed
660
661
662
#if defined (NATIVE_PRECISION)
        partial_energies[get_local_id(0)] += native_divide (
                                                             dockpars_coeff_elec * atom_charges_const[atom1_id] * atom_charges_const[atom2_id],
663
                                                             atomic_distance * (-8.5525f + native_divide(86.9525f,(1.0f + 7.7839f*native_exp(-0.3154f*atomic_distance))))
Leonardo Solis's avatar
Leonardo Solis committed
664
665
666
667
                                                             );
#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],
668
                                                             atomic_distance * (-8.5525f + half_divide(86.9525f,(1.0f + 7.7839f*half_exp(-0.3154f*atomic_distance))))
Leonardo Solis's avatar
Leonardo Solis committed
669
670
                                                             );
#else	// Full precision
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
				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]/
689
			                                       (atomic_distance*(-8.5525f + 86.9525f/(1.0f + 7.7839f*exp(-0.3154f*atomic_distance))));
Leonardo Solis's avatar
Leonardo Solis committed
690
#endif
Leonardo Solis's avatar
Leonardo Solis committed
691
692
693
694
695
696
697
698
*/


#if defined (NATIVE_PRECISION)
        partial_energies[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))))
                                                             );
699
700
701
702
703
704
705
706
707

			#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


Leonardo Solis's avatar
Leonardo Solis committed
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
#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],
                                                             atomic_distance * (DIEL_A + half_divide(DIEL_B,(1.0f + DIEL_K*half_exp(-DIEL_B_TIMES_H*atomic_distance))))
                                                             );
#else	// Full precision
	partial_energies[get_local_id(0)] += dockpars_coeff_elec*atom_charges_const[atom1_id]*atom_charges_const[atom2_id]/
			                                       (atomic_distance*(DIEL_A + DIEL_B/(1.0f + DIEL_K*exp(-DIEL_B_TIMES_H*atomic_distance))));
#endif







Leonardo Solis's avatar
Leonardo Solis committed
724

Leonardo Solis's avatar
Leonardo Solis committed
725
			// Calculating desolvation term
Leonardo Solis's avatar
Leonardo Solis committed
726
727
#if defined (NATIVE_PRECISION)
			partial_energies[get_local_id(0)] += ((dspars_S_const[atom1_typeid] +
728
729
730
731
							       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));
732
733
734
735
736
737
738
739
740

			#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

Leonardo Solis's avatar
Leonardo Solis committed
741
742
#elif defined (HALF_PRECISION)
			partial_energies[get_local_id(0)] += ((dspars_S_const[atom1_typeid] +
743
744
745
746
							       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*half_exp(-atomic_distance*half_divide(atomic_distance,25.92f));
Leonardo Solis's avatar
Leonardo Solis committed
747
748
#else	// Full precision
			partial_energies[get_local_id(0)] += ((dspars_S_const[atom1_typeid] +
749
750
751
752
							       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*exp(-atomic_distance*atomic_distance/25.92f);
Leonardo Solis's avatar
Leonardo Solis committed
753
754
755
#endif

		}
756
757
	} // End contributor_counter for-loop (INTRAMOLECULAR ENERGY)

Leonardo Solis's avatar
Leonardo Solis committed
758
759
760
761
762
763
	barrier(CLK_LOCAL_MEM_FENCE);

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

764
765
766
		for (uint contributor_counter=1;
		          contributor_counter<NUM_OF_THREADS_PER_BLOCK;
		          contributor_counter++)
Leonardo Solis's avatar
Leonardo Solis committed
767
768
769
770
		{
			*energy += partial_energies[contributor_counter];
		}
	}
771

772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
	#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

Leonardo Solis's avatar
Leonardo Solis committed
792
793
794
795
796
797
798
}

#include "kernel1.cl"
#include "kernel2.cl"
#include "auxiliary_genetic.cl"
#include "kernel4.cl"
#include "kernel3.cl"
799
#include "auxiliary_gradient.cl"
800
#include "calcgradient.cl"
801
#include "kernel_gradient.cl"