calcenergy.cl 23.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.

*/

24
25
26
//#define DEBUG_ENERGY_KERNEL1
//#define DEBUG_ENERGY_KERNEL4
//#define DEBUG_ENERGY_KERNEL3
27

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
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,
38
39
40
								    		// 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
41
42
43
44
45
46
47
		 __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
48

49
50
51
52
                    // 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.
53
54
55
56
57
58
59
60
61
		    	__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,

62
			#if defined (DEBUG_ENERGY_KERNEL1) || defined (DEBUG_ENERGY_KERNEL4) || defined (DEBUG_ENERGY_KERNEL3)
63
64
65
66
			__local float* partial_interE,
			__local float* partial_intraE,
			#endif

67
68
69
	             __constant float* atom_charges_const,
                     __constant char*  atom_types_const,
                     __constant char*  intraE_contributors_const,
lvs's avatar
lvs committed
70
71
72
	                  	float  dockpars_smooth,
	       	     __constant float* reqm,
	       	     __constant float* reqm_hbond,
73
74
75
76
77
78
79
80
81
82
                     __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,
83
                     __constant float* ref_orientation_quats_const
Leonardo Solis's avatar
Leonardo Solis committed
84
85
86
87
88
89
90
91
92
)

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

93
	#if defined (DEBUG_ENERGY_KERNEL1) || defined (DEBUG_ENERGY_KERNEL4) || defined (DEBUG_ENERGY_KERNEL3)
94
95
96
97
	partial_interE[get_local_id(0)] = 0.0f;
	partial_intraE[get_local_id(0)] = 0.0f;
	#endif

98
	uchar g1 = dockpars_gridsize_x;
99
100
	uint  g2 = dockpars_gridsize_x_times_y /*dockpars_gridsize_x * dockpars_gridsize_y*/;
  	uint  g3 = dockpars_gridsize_x_times_y_times_z /*dockpars_gridsize_x * dockpars_gridsize_y * dockpars_gridsize_z*/;
101

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

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

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

118
			if ((rotation_list_element & RLIST_FIRSTROT_MASK) != 0)	// If first rotation of this atom
Leonardo Solis's avatar
Leonardo Solis committed
119
120
121
122
123
124
125
126
127
128
129
130
			{
				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];
			}

131
132
133
134
135
136
137
			// 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
138
			{
139
				// Rotational genes in the Shoemake space expressed in radians
140
				float u1 = genotype[3];
Leonardo Solis's avatar
Leonardo Solis committed
141
142
				float u2 = genotype[4];
				float u3 = genotype[5];
143
144
145
146
147
148

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

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

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

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

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

				// Transforming torsion angles into quaternions
175
176
177
178
179
180
				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
181
182
			}

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

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

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

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

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

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

331
			#if defined (DEBUG_ENERGY_KERNEL1) || defined (DEBUG_ENERGY_KERNEL4) || defined (DEBUG_ENERGY_KERNEL3)
332
333
334
			partial_interE[get_local_id(0)] += TRILININTERPOL(cube, weights);
			#endif

Leonardo Solis's avatar
Leonardo Solis committed
335
			// Capturing electrostatic values
336
			atom_typeid = dockpars_num_of_atypes;
Leonardo Solis's avatar
Leonardo Solis committed
337

338
339
340
341
342
343
344
345
346
			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);
347
348

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

351
			#if defined (DEBUG_ENERGY_KERNEL1) || defined (DEBUG_ENERGY_KERNEL4) || defined (DEBUG_ENERGY_KERNEL3)
352
353
354
			partial_interE[get_local_id(0)] += q * TRILININTERPOL(cube, weights);
			#endif

Leonardo Solis's avatar
Leonardo Solis committed
355
			// Capturing desolvation values
356
			atom_typeid = dockpars_num_of_atypes+1;
Leonardo Solis's avatar
Leonardo Solis committed
357

358
359
360
361
362
363
364
365
366
			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
367

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

371
			#if defined (DEBUG_ENERGY_KERNEL1) || defined (DEBUG_ENERGY_KERNEL4) || defined (DEBUG_ENERGY_KERNEL3)
372
373
			partial_interE[get_local_id(0)] += fabs(q) * TRILININTERPOL(cube, weights);
			#endif
Leonardo Solis's avatar
Leonardo Solis committed
374
375
		}

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

378

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

	// ================================================
Leonardo Solis's avatar
Leonardo Solis committed
406
	// CALCULATING INTRAMOLECULAR ENERGY
Leonardo Solis's avatar
Leonardo Solis committed
407
	// ================================================
408
409
410
	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
411
	{
Leonardo Solis's avatar
Leonardo Solis committed
412
		// Getting atom IDs
413
414
		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
415

416
417
418
419
420
		// 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
421

Leonardo Solis's avatar
Leonardo Solis committed
422
		// Calculating atomic_distance
423
		float atomic_distance = native_sqrt(subx*subx + suby*suby + subz*subz)*dockpars_grid_spacing;
Leonardo Solis's avatar
Leonardo Solis committed
424

Leonardo Solis's avatar
Leonardo Solis committed
425
		// Calculating energy contributions
lvs's avatar
lvs committed
426
		if (atomic_distance < 8.0f)
Leonardo Solis's avatar
Leonardo Solis committed
427
		{
428
429
430
			// 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
431

lvs's avatar
lvs committed
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
			// 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
			{
				opt_distance = reqm_hbond [atom1_typeid] + reqm_hbond [atom2_typeid];
			}
			else	//van der Waals
			{
				opt_distance = 0.5f*(reqm [atom1_typeid] + reqm [atom2_typeid]);
			}

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

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

				if (intraE_contributors_const[3*contributor_counter+2] == 1)	//H-bond
				{
					printf("%-5s %u %u %f %f %f %f %f %f\n", "hbond", atom1_typeid, atom2_typeid, reqm_hbond [atom1_typeid], reqm_hbond [atom2_typeid], opt_distance, delta_distance, atomic_distance, smoothed_distance);

				}
				else	//van der Waals
				{
					printf("%-5s %u %u %f %f %f %f %f %f\n", "vdw", atom1_typeid, atom2_typeid, reqm [atom1_typeid], reqm [atom2_typeid], opt_distance, delta_distance, atomic_distance, smoothed_distance);	
				}
			}
*/

Leonardo Solis's avatar
Leonardo Solis committed
478
			// Calculating van der Waals / hydrogen bond term
lvs's avatar
lvs committed
479
			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));
480

481
			#if defined (DEBUG_ENERGY_KERNEL1) || defined (DEBUG_ENERGY_KERNEL4) || defined (DEBUG_ENERGY_KERNEL3)
lvs's avatar
lvs committed
482
			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));
483
484
485
			#endif

			if (intraE_contributors_const[3*contributor_counter+2] == 1) {	//H-bond
lvs's avatar
lvs committed
486
				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));
487
				#if defined (DEBUG_ENERGY_KERNEL1) || defined (DEBUG_ENERGY_KERNEL4) || defined (DEBUG_ENERGY_KERNEL3)
lvs's avatar
lvs committed
488
				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));
489
				#endif
490
491
			}
			else {	//van der Waals
lvs's avatar
lvs committed
492
				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));
493

494
				#if defined (DEBUG_ENERGY_KERNEL1) || defined (DEBUG_ENERGY_KERNEL4) || defined (DEBUG_ENERGY_KERNEL3)
lvs's avatar
lvs committed
495
				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));
496
				#endif
497
			}
lvs's avatar
lvs committed
498

Leonardo Solis's avatar
Leonardo Solis committed
499
			// Calculating electrostatic term
500
        		partial_energies[get_local_id(0)] += native_divide (
Leonardo Solis's avatar
Leonardo Solis committed
501
502
503
                                                             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))))
                                                             );
504

505
			#if defined (DEBUG_ENERGY_KERNEL1) || defined (DEBUG_ENERGY_KERNEL4) || defined (DEBUG_ENERGY_KERNEL3)
506
507
508
509
510
511
			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
512
			// Calculating desolvation term
513
			// 1/25.92 = 0.038580246913580245
Leonardo Solis's avatar
Leonardo Solis committed
514
			partial_energies[get_local_id(0)] += ((dspars_S_const[atom1_typeid] +
515
516
517
							       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]) *
518
					                       dockpars_coeff_desolv*native_exp(-0.03858025f*native_powr(atomic_distance, 2));
519

520
			#if defined (DEBUG_ENERGY_KERNEL1) || defined (DEBUG_ENERGY_KERNEL4) || defined (DEBUG_ENERGY_KERNEL3)
521
522
523
524
			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]) *
525
					                       dockpars_coeff_desolv*native_exp(-0.03858025f*native_powr(atomic_distance, 2));
526
527
			#endif

528

Leonardo Solis's avatar
Leonardo Solis committed
529
530

		}
531
532
	} // End contributor_counter for-loop (INTRAMOLECULAR ENERGY)

Leonardo Solis's avatar
Leonardo Solis committed
533
534
535
536
537
538
	barrier(CLK_LOCAL_MEM_FENCE);

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

539
540
541
		for (uint contributor_counter=1;
		          contributor_counter<NUM_OF_THREADS_PER_BLOCK;
		          contributor_counter++)
Leonardo Solis's avatar
Leonardo Solis committed
542
543
544
545
		{
			*energy += partial_energies[contributor_counter];
		}
	}
546

547
	#if defined (DEBUG_ENERGY_KERNEL1) || defined (DEBUG_ENERGY_KERNEL4) || defined (DEBUG_ENERGY_KERNEL3)
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
	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
567
568
569
570
571
572
573
}

#include "kernel1.cl"
#include "kernel2.cl"
#include "auxiliary_genetic.cl"
#include "kernel4.cl"
#include "kernel3.cl"
574
//#include "auxiliary_gradient.cl"
575
#include "calcgradient.cl"
576
#include "kernel_gradient.cl"