calcenergy.cl 26.9 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
#include "calcenergy_basic.h"
Leonardo Solis's avatar
Leonardo Solis committed
27
28
29

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

30
31
32
33
34
35
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,
36
37
38
								    		// 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
39
40
41
42
43
44
45
		 __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
46

47
48
49
50
                    // 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.
51
52
53
54
55
56
57
58
59
		    	__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
60
			#if defined (DEBUG_ENERGY_KERNEL)
61
62
63
64
			__local float* partial_interE,
			__local float* partial_intraE,
			#endif

65
66
67
	             __constant float* atom_charges_const,
                     __constant char*  atom_types_const,
                     __constant char*  intraE_contributors_const,
lvs's avatar
lvs committed
68
69
70
#if 0
 				bool   debug,
#endif
lvs's avatar
lvs committed
71
72
73
	                  	float  dockpars_smooth,
	       	     __constant float* reqm,
	       	     __constant float* reqm_hbond,
lvs's avatar
lvs committed
74
75
	     	     __constant uint*  atom1_types_reqm,
	     	     __constant uint*  atom2_types_reqm,
76
77
78
79
80
81
82
83
84
85
                     __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,
86
                     __constant float* ref_orientation_quats_const
Leonardo Solis's avatar
Leonardo Solis committed
87
88
89
90
91
92
93
94
95
)

//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
96
	#if defined (DEBUG_ENERGY_KERNEL)
97
98
99
100
	partial_interE[get_local_id(0)] = 0.0f;
	partial_intraE[get_local_id(0)] = 0.0f;
	#endif

101
	uchar g1 = dockpars_gridsize_x;
lvs's avatar
lvs committed
102
103
	uint  g2 = dockpars_gridsize_x_times_y;
  	uint  g3 = dockpars_gridsize_x_times_y_times_z;
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
				// Rotational genes in the Shoemake space expressed in radians
143
				float u1 = genotype[3];
Leonardo Solis's avatar
Leonardo Solis committed
144
145
				float u2 = genotype[4];
				float u3 = genotype[5];
146
147
148
149
150
151

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

Leonardo Solis's avatar
Leonardo Solis committed
153
154
155
156
				rotation_movingvec[0] = genotype[0];
				rotation_movingvec[1] = genotype[1];
				rotation_movingvec[2] = genotype[2];
			}
157
			else	// If rotating around rotatable bond
Leonardo Solis's avatar
Leonardo Solis committed
158
			{
159
				uint rotbond_id = (rotation_list_element & RLIST_RBONDID_MASK) >> RLIST_RBONDID_SHIFT;
Leonardo Solis's avatar
Leonardo Solis committed
160

161
				float rotation_unitvec[3];
Leonardo Solis's avatar
Leonardo Solis committed
162
163
164
				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];
165
				float rotation_angle = genotype[6+rotbond_id]*DEG_TO_RAD;
Leonardo Solis's avatar
Leonardo Solis committed
166
167
168
169
170

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

171
172
				// Performing additionally the first movement which 
				// is needed only if rotating around rotatable bond
Leonardo Solis's avatar
Leonardo Solis committed
173
174
175
				atom_to_rotate[0] -= rotation_movingvec[0];
				atom_to_rotate[1] -= rotation_movingvec[1];
				atom_to_rotate[2] -= rotation_movingvec[2];
176
177

				// Transforming torsion angles into quaternions
lvs's avatar
lvs committed
178
				rotation_angle  = rotation_angle * 0.5f;
179
180
181
182
183
				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
184
185
			}

186
187
188
189
			// 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
190
			{
191
192
				// Calculating quatrot_left*ref_orientation_quats_const,
				// which means that reference orientation rotation is the first
Leonardo Solis's avatar
Leonardo Solis committed
193
194
195
196
197
198
				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)]-
199
200
201
						 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
202
				quatrot_left_x = quatrot_temp_q*ref_orientation_quats_const[4*(*run_id)+1]+
203
204
205
						 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
206
				quatrot_left_y = quatrot_temp_q*ref_orientation_quats_const[4*(*run_id)+2]+
207
208
209
						 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
210
				quatrot_left_z = quatrot_temp_q*ref_orientation_quats_const[4*(*run_id)+3]+
211
212
213
						 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
214
215
216
			}

			quatrot_temp_q = 0 -
217
218
219
					 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
220
			quatrot_temp_x = quatrot_left_q*atom_to_rotate [0] +
221
222
					 quatrot_left_y*atom_to_rotate [2] -
					 quatrot_left_z*atom_to_rotate [1];
Leonardo Solis's avatar
Leonardo Solis committed
223
			quatrot_temp_y = quatrot_left_q*atom_to_rotate [1] -
224
225
					 quatrot_left_x*atom_to_rotate [2] +
					 quatrot_left_z*atom_to_rotate [0];
Leonardo Solis's avatar
Leonardo Solis committed
226
			quatrot_temp_z = quatrot_left_q*atom_to_rotate [2] +
227
228
					 quatrot_left_x*atom_to_rotate [1] -
					 quatrot_left_y*atom_to_rotate [0];
Leonardo Solis's avatar
Leonardo Solis committed
229
230

			atom_to_rotate [0] = 0 -
231
232
233
234
					  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
235
			atom_to_rotate [1] = 0 -
236
237
238
239
					  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
240
			atom_to_rotate [2] = 0 -
241
242
243
244
					  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
245

246
			// Performing final movement and storing values
Leonardo Solis's avatar
Leonardo Solis committed
247
248
249
250
251
252
253
254
255
256
257
			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
258
	// CALCULATING INTERMOLECULAR ENERGY
Leonardo Solis's avatar
Leonardo Solis committed
259
	// ================================================
260
261
262
	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
263
	{
264
265
266
267
268
		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
269
270

		if ((x < 0) || (y < 0) || (z < 0) || (x >= dockpars_gridsize_x-1)
271
272
				                  || (y >= dockpars_gridsize_y-1)
						  || (z >= dockpars_gridsize_z-1)){
Leonardo Solis's avatar
Leonardo Solis committed
273
			partial_energies[get_local_id(0)] += 16777216.0f; //100000.0f;
274
	
lvs's avatar
lvs committed
275
			#if defined (DEBUG_ENERGY_KERNEL)
276
277
			partial_interE[get_local_id(0)] += 16777216.0f;
			#endif
Leonardo Solis's avatar
Leonardo Solis committed
278
279
280
		}
		else
		{
Leonardo Solis's avatar
Leonardo Solis committed
281
			// Getting coordinates
282
283
284
285
286
287
288
289
290
291
			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
292
			// Calculating interpolation weights
293
			float weights[2][2][2];
Leonardo Solis's avatar
Leonardo Solis committed
294
295
296
297
298
299
300
301
302
			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
303
			// Capturing affinity values
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
			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
330

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

lvs's avatar
lvs committed
334
			#if defined (DEBUG_ENERGY_KERNEL)
335
336
337
			partial_interE[get_local_id(0)] += TRILININTERPOL(cube, weights);
			#endif

Leonardo Solis's avatar
Leonardo Solis committed
338
			// Capturing electrostatic values
339
			atom_typeid = dockpars_num_of_atypes;
Leonardo Solis's avatar
Leonardo Solis committed
340

341
342
343
344
345
346
347
348
349
			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);
350
351

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

lvs's avatar
lvs committed
354
			#if defined (DEBUG_ENERGY_KERNEL)
355
356
357
			partial_interE[get_local_id(0)] += q * TRILININTERPOL(cube, weights);
			#endif

Leonardo Solis's avatar
Leonardo Solis committed
358
			// Capturing desolvation values
359
			atom_typeid = dockpars_num_of_atypes+1;
Leonardo Solis's avatar
Leonardo Solis committed
360

361
362
363
364
365
366
367
368
369
			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
370

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

lvs's avatar
lvs committed
374
			#if defined (DEBUG_ENERGY_KERNEL)
375
376
			partial_interE[get_local_id(0)] += fabs(q) * TRILININTERPOL(cube, weights);
			#endif
Leonardo Solis's avatar
Leonardo Solis committed
377
378
		}

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

381

lvs's avatar
lvs committed
382
	#if defined (DEBUG_ENERGY_KERNEL)
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
	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
403
404
	// In paper: intermolecular and internal energy calculation
	// are independent from each other, -> NO BARRIER NEEDED
405
  	// but require different operations,
Leonardo Solis's avatar
Leonardo Solis committed
406
407
408
	// thus, they can be executed only sequentially on the GPU.

	// ================================================
Leonardo Solis's avatar
Leonardo Solis committed
409
	// CALCULATING INTRAMOLECULAR ENERGY
Leonardo Solis's avatar
Leonardo Solis committed
410
	// ================================================
411
412
413
	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
414
415
416
417
418
419
420

#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
421
	{
lvs's avatar
lvs committed
422
423
424
425
426
427
428
429
430
431
#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
432
		// Getting atom IDs
433
434
		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
435

436
437
438
439
440
		// 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
441

Leonardo Solis's avatar
Leonardo Solis committed
442
		// Calculating atomic_distance
443
		float atomic_distance = native_sqrt(subx*subx + suby*suby + subz*subz)*dockpars_grid_spacing;
Leonardo Solis's avatar
Leonardo Solis committed
444

Leonardo Solis's avatar
Leonardo Solis committed
445
		// Calculating energy contributions
lvs's avatar
lvs committed
446
/*
lvs's avatar
lvs committed
447
		if (atomic_distance < 8.0f)
Leonardo Solis's avatar
Leonardo Solis committed
448
		{
lvs's avatar
lvs committed
449
*/
450
451
452
			// 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
453

lvs's avatar
lvs committed
454
455
456
			uint atom1_type_vdw_hb = atom1_types_reqm [atom1_typeid];
	     	        uint atom2_type_vdw_hb = atom2_types_reqm [atom2_typeid];

lvs's avatar
lvs committed
457
458
459
460
461
462
463
464
465
			// Getting optimum pair distance (opt_distance) from reqm and reqm_hbond
			// reqm: equilibrium internuclear separation 
			//       (sum of the vdW radii of two like atoms (A)) in the case of vdW
			// reqm_hbond: equilibrium internuclear separation
			//  	 (sum of the vdW radii of two like atoms (A)) in the case of hbond 
			float opt_distance;

			if (intraE_contributors_const[3*contributor_counter+2] == 1)	//H-bond
			{
lvs's avatar
lvs committed
466
#if 0
lvs's avatar
lvs committed
467
				opt_distance = reqm_hbond [atom1_typeid] + reqm_hbond [atom2_typeid];
lvs's avatar
lvs committed
468
469
#endif
				opt_distance = reqm_hbond [atom1_type_vdw_hb] + reqm_hbond [atom2_type_vdw_hb];
lvs's avatar
lvs committed
470
471
472
			}
			else	//van der Waals
			{
lvs's avatar
lvs committed
473
#if 0
lvs's avatar
lvs committed
474
				opt_distance = 0.5f*(reqm [atom1_typeid] + reqm [atom2_typeid]);
lvs's avatar
lvs committed
475
476
#endif
				opt_distance = 0.5f*(reqm [atom1_type_vdw_hb] + reqm [atom2_type_vdw_hb]);
lvs's avatar
lvs committed
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
			}

			// Getting smoothed distance
			// smoothed_distance = function(atomic_distance, opt_distance)
			float smoothed_distance;
			float delta_distance = 0.5f*dockpars_smooth;

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

Leonardo Solis's avatar
Leonardo Solis committed
494
			// Calculating van der Waals / hydrogen bond term
lvs's avatar
lvs committed
495
			partial_energies[get_local_id(0)] += native_divide(VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance/*atomic_distance*/,12));
lvs's avatar
lvs committed
496
497
498
499
500
#if 0
smoothed_intraE = native_divide(VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance,12));
raw_intraE_vdw_hb      = native_divide(VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance  ,12)); 
#endif
			#if defined (DEBUG_ENERGY_KERNEL)
lvs's avatar
lvs committed
501
			partial_intraE[get_local_id(0)] += native_divide(VWpars_AC_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance/*atomic_distance*/,12));
502
503
504
			#endif

			if (intraE_contributors_const[3*contributor_counter+2] == 1) {	//H-bond
lvs's avatar
lvs committed
505
				partial_energies[get_local_id(0)] -= native_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance/*atomic_distance*/,10));
lvs's avatar
lvs committed
506
507
508
509
510
#if 0
smoothed_intraE -= native_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance,10));
raw_intraE_vdw_hb 	-= native_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance  ,10));
#endif
				#if defined (DEBUG_ENERGY_KERNEL)
lvs's avatar
lvs committed
511
				partial_intraE[get_local_id(0)] -= native_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance/*atomic_distance*/,10));
512
				#endif
513
514
			}
			else {	//van der Waals
lvs's avatar
lvs committed
515
				partial_energies[get_local_id(0)] -= native_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance/*atomic_distance*/,6));
lvs's avatar
lvs committed
516
517
518
519
520
#if 0
smoothed_intraE -= native_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance,6));
raw_intraE_vdw_hb      -= native_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(atomic_distance  ,6));
#endif
				#if defined (DEBUG_ENERGY_KERNEL)
lvs's avatar
lvs committed
521
				partial_intraE[get_local_id(0)] -= native_divide(VWpars_BD_const[atom1_typeid * dockpars_num_of_atypes+atom2_typeid],native_powr(smoothed_distance/*atomic_distance*/,6));
522
				#endif
523
			}
lvs's avatar
lvs committed
524

Leonardo Solis's avatar
Leonardo Solis committed
525
			// Calculating electrostatic term
526
        		partial_energies[get_local_id(0)] += native_divide (
Leonardo Solis's avatar
Leonardo Solis committed
527
528
529
                                                             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))))
                                                             );
lvs's avatar
lvs committed
530
531
532
533
534
535
536
537
538
539
540
541
542
#if 0
smoothed_intraE += 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))))
                             );

raw_intraE_el 	= 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

			#if defined (DEBUG_ENERGY_KERNEL)
543
544
545
546
547
548
			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
549
			// Calculating desolvation term
550
			// 1/25.92 = 0.038580246913580245
Leonardo Solis's avatar
Leonardo Solis committed
551
			partial_energies[get_local_id(0)] += ((dspars_S_const[atom1_typeid] +
552
553
554
							       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]) *
555
					                       dockpars_coeff_desolv*native_exp(-0.03858025f*native_powr(atomic_distance, 2));
556

lvs's avatar
lvs committed
557
558
559
560
561
562
563
564
565
566
567
568
569
570
#if 0
smoothed_intraE += ((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(-0.03858025f*native_powr(atomic_distance, 2));

raw_intraE_sol = ((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(-0.03858025f*native_powr(atomic_distance, 2));
#endif
			#if defined (DEBUG_ENERGY_KERNEL)
571
572
573
574
			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]) *
575
					                       dockpars_coeff_desolv*native_exp(-0.03858025f*native_powr(atomic_distance, 2));
576
577
			#endif

578

Leonardo Solis's avatar
Leonardo Solis committed
579

lvs's avatar
lvs committed
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628




#if 0
			raw_intraE = raw_intraE_vdw_hb + raw_intraE_el + raw_intraE_sol;

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

				if (intraE_contributors_const[3*contributor_counter+2] == 1)	//H-bond
				{
					// diogos table
//					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, reqm_hbond [atom1_type_vdw_hb], reqm_hbond [atom2_type_vdw_hb], opt_distance, atomic_distance, smoothed_distance, smoothed_intraE, raw_intraE);

					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, reqm_hbond [atom1_type_vdw_hb], reqm_hbond [atom2_type_vdw_hb], opt_distance, atomic_distance, smoothed_distance, smoothed_intraE, raw_intraE, raw_intraE_vdw_hb, raw_intraE_el, raw_intraE_sol);




				}
				else	//van der Waals
				{
					// diogos table
//					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, reqm [atom1_type_vdw_hb], reqm [atom2_type_vdw_hb], opt_distance, atomic_distance, smoothed_distance, smoothed_intraE, raw_intraE);

					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, reqm [atom1_type_vdw_hb], reqm [atom2_type_vdw_hb], opt_distance, atomic_distance, smoothed_distance, smoothed_intraE, raw_intraE, raw_intraE_vdw_hb, raw_intraE_el, raw_intraE_sol);


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









/*
		} // if (atomic_distance < 8.0f)
*/


629
630
	} // End contributor_counter for-loop (INTRAMOLECULAR ENERGY)

lvs's avatar
lvs committed
631
632
633
634
#if 0
} // if (get_local_id (0) == 0) {
#endif

Leonardo Solis's avatar
Leonardo Solis committed
635
636
637
638
639
640
	barrier(CLK_LOCAL_MEM_FENCE);

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

641
642
643
		for (uint contributor_counter=1;
		          contributor_counter<NUM_OF_THREADS_PER_BLOCK;
		          contributor_counter++)
Leonardo Solis's avatar
Leonardo Solis committed
644
645
646
647
		{
			*energy += partial_energies[contributor_counter];
		}
	}
648

649
650
	barrier(CLK_LOCAL_MEM_FENCE);

lvs's avatar
lvs committed
651
	#if defined (DEBUG_ENERGY_KERNEL)
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
	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
669
670
671
672
673
674
675
}

#include "kernel1.cl"
#include "kernel2.cl"
#include "auxiliary_genetic.cl"
#include "kernel4.cl"
#include "kernel3.cl"
676
#include "calcgradient.cl"
677
#include "kernel_gradient.cl"